HyperLib 1997 Winter - Disc 1/HYPERLIB-1997-Winter-CD1.ISO.7z/HYPERLIB-1997-Winter-CD1.ISO/オンラインウェア/PRG/PL 2.0 SupplementDoc Folder.sit/PL 2.0 SupplementDoc Folder/Documentation/Chapter 12. Numbers

Text File | 1995-03-27 | 132KB | 3,162 lines

Common Lisp the Language, 2nd Edition
-------------------------------------------------------------------------------

12. Numbers

Common Lisp provides several different representations for numbers. These
representations may be divided into four categories: integers, ratios,
floating-point numbers, and complex numbers. Many numeric functions will accept
any kind of number; they are generic. Other functions accept only certain kinds
of numbers.

[change_begin]
Note that this remark, predating the design of the Common Lisp Object System,
uses the term ``generic'' in a generic sense and not necessarily in the
technical sense used by CLOS (see chapter 2).
[change_end]

In general, numbers in Common Lisp are not true objects; eq cannot be counted
upon to operate on them reliably. In particular, it is possible that the
expression

(let ((x z) (y z)) (eq x y))

may be false rather than true if the value of z is a number.

-------------------------------------------------------------------------------
Rationale: This odd breakdown of eq in the case of numbers allows the
implementor enough design freedom to produce exceptionally efficient numerical
code on conventional architectures. MacLisp requires this freedom, for example,
in order to produce compiled numerical code equal in speed to Fortran. Common
Lisp makes this same restriction, if not for this freedom, then at least for
the sake of compatibility.
-------------------------------------------------------------------------------

If two objects are to be compared for ``identity,'' but either might be a
number, then the predicate eql is probably appropriate; if both objects are
known to be numbers, then = may be preferable.

-------------------------------------------------------------------------------

   *  Precision, Contagion, and Coercion
   *  Predicates on Numbers
   *  Comparisons on Numbers
   *  Arithmetic Operations
   *  Irrational and Transcendental Functions
        o  Exponential and Logarithmic Functions
        o  Trigonometric and Related Functions
        o  Branch Cuts, Principal Values, and Boundary Conditions in the
          Complex Plane
   *  Type Conversions and Component Extractions on Numbers
   *  Logical Operations on Numbers
   *  Byte Manipulation Functions
   *  Random Numbers
   *  Implementation Parameters

-------------------------------------------------------------------------------

12.1. Precision, Contagion, and Coercion

In general, computations with floating-point numbers are only approximate. The
precision of a floating-point number is not necessarily correlated at all with
the accuracy of that number. For instance, 3.142857142857142857 is a more
precise approximation to than 3.14159, but the latter is more accurate. The
precision refers to the number of bits retained in the representation. When an
operation combines a short floating-point number with a long one, the result
will be a long floating-point number. This rule is made to ensure that as much
accuracy as possible is preserved; however, it is by no means a guarantee.
Common Lisp numerical routines do assume, however, that the accuracy of an
argument does not exceed its precision. Therefore when two small floating-point
numbers are combined, the result will always be a small floating-point number.
This assumption can be overridden by first explicitly converting a small
floating-point number to a larger representation. (Common Lisp never converts
automatically from a larger size to a smaller one.)

Rational computations cannot overflow in the usual sense (though of course
there may not be enough storage to represent one), as integers and ratios may
in principle be of any magnitude. Floating-point computations may get exponent
overflow or underflow; this is an error.

[change_begin]
X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to address certain problems
relating to floating-point overflow and underflow, but certain parts of the
proposed solution were not adopted, namely to add the macro
without-floating-underflow-traps to the language and to require certain
behavior of floating-point overflow and underflow. The committee agreed that
this area of the language requires more discussion before a solution is
standardized.

For the record, the proposal that was considered and rejected (for the nonce)
introduced a macro without-floating-underflow-traps that would execute its body
in such a way that, within its dynamic extent, a floating-point underflow must
not signal an error but instead must produce either a denormalized number or
zero as the result. The rejected proposal also specified the following
treatment of overflow and underflow:

   *  A floating-point computation that overflows should signal an error of
     type floating-point-overflow.

   *  Unless the dynamic extent of a use of without-floating-underflow-traps, a
     floating-point computation that underflows should signal an error of type
     floating-point-underflow. A result that can be represented only in
     denormalized form must be considered an underflow in implementations that
     support denormalized floating-point numbers.

These points refer to conditions floating-point-overflow and
floating-point-underflow that were approved by X3J13 and are described in
section 29.5.
[change_end]

When rational and floating-point numbers are compared or combined by a
numerical function, the rule of floating-point contagion is followed: when a
rational meets a floating-point number, the rational is first converted to a
floating-point number of the same format. For functions such as + that take
more than two arguments, it may be that part of the operation is carried out
exactly using rationals and then the rest is done using floating-point
arithmetic.

[change_begin]
X3J13 voted in January 1989 (CONTAGION-ON-NUMERICAL-COMPARISONS)   to apply the
rule of floating-point contagion stated above to the case of combining rational
and floating-point numbers. For comparing, the following rule is to be used
instead: When a rational number and a floating-point number are to be compared
by a numerical function, in effect the floating-point number is first converted
to a rational number as if by the function rational, and then an exact
comparison of two rational numbers is performed. It is of course valid to use a
more efficient implementation than actually calling the function rational, as
long as the result of the comparison is the same. In the case of complex
numbers, the real and imaginary parts are handled separately.

-------------------------------------------------------------------------------
Rationale: In general, accuracy cannot be preserved in combining operations,
but it can be preserved in comparisons, and preserving it makes that part of
Common Lisp algebraically a bit more tractable. In particular, this change
prevents the breakdown of transitivity. Let a be the result of (/ 10.0
single-float-epsilon), and let j be the result of (floor a). (Note that (= a (+
a 1.0)) is true, by the definition of single-float-epsilon.) Under the old
rules, all of (<= a j), (< j (+ j 1)), and (<= (+ j 1) a) would be true;
transitivity would then imply that (< a a) ought to be true, but of course it
is false, and therefore transitivity fails. Under the new rule, however, (<= (+
j 1) a) is false.
-------------------------------------------------------------------------------

[change_end]

For functions that are mathematically associative (and possibly commutative), a
Common Lisp implementation may process the arguments in any manner consistent
with associative (and possibly commutative) rearrangement. This does not affect
the order in which the argument forms are evaluated, of course; that order is
always left to right, as in all Common Lisp function calls. What is left loose
is the order in which the argument values are processed. The point of all this
is that implementations may differ in which automatic coercions are applied
because of differing orders of argument processing. As an example, consider
this expression:

(+ 1/3 2/3 1.0D0 1.0 1.0E-15)

One implementation might process the arguments from left to right, first adding
1/3 and 2/3 to get 1, then converting that to a double-precision floating-point
number for combination with 1.0D0, then successively converting and adding 1.0
and 1.0E-15. Another implementation might process the arguments from right to
left, first performing a single-precision floating-point addition of 1.0 and
1.0E-15 (and probably losing some accuracy in the process!), then converting
the sum to double precision and adding 1.0D0, then converting 2/3 to
double-precision floating-point and adding it, and then converting 1/3 and
adding that. A third implementation might first scan all the arguments, process
all the rationals first to keep that part of the computation exact, then find
an argument of the largest floating-point format among all the arguments and
add that, and then add in all other arguments, converting each in turn (all in
a perhaps misguided attempt to make the computation as accurate as possible).
In any case, all three strategies are legitimate. The user can of course
control the order of processing explicitly by writing several calls; for
example:

(+ (+ 1/3 2/3) (+ 1.0D0 1.0E-15) 1.0)

The user can also control all coercions simply by writing calls to coercion
functions explicitly.

In general, then, the type of the result of a numerical function is a
floating-point number of the largest format among all the floating-point
arguments to the function; but if the arguments are all rational, then the
result is rational (except for functions that can produce mathematically
irrational results, in which case a single-format floating-point number may
result).

There is a separate rule of complex contagion. As a rule, complex numbers never
result from a numerical function unless one or more of the arguments is
complex. (Exceptions to this rule occur among the irrational and transcendental
functions, specifically expt, log, sqrt, asin, acos, acosh, and atanh; see
section 12.5.) When a non-complex number meets a complex number, the
non-complex number is in effect first converted to a complex number by
providing an imaginary part of zero.

If any computation produces a result that is a ratio of two integers such that
the denominator evenly divides the numerator, then the result is immediately
converted to the equivalent integer. This is called the rule of rational
canonicalization.

If the result of any computation would be a complex rational with a zero
imaginary part, the result is immediately converted to a non-complex rational
number by taking the real part. This is called the rule of complex
canonicalization. Note that this rule does not apply to complex numbers whose
components are floating-point numbers. Whereas #C(5 0) and 5 are not distinct
values in Common Lisp (they are always eql), #C(5.0 0.0) and 5.0 are always
distinct values in Common Lisp (they are never eql, although they are equalp).

-------------------------------------------------------------------------------

12.2. Predicates on Numbers

Each of the following functions tests a single number for a specific property.
Each function requires that its argument be a number; to call one with a
non-number is an error.

[Function]
zerop number

This predicate is true if number is zero (the integer zero, a floating-point
zero, or a complex zero), and is false otherwise. Regardless of whether an
implementation provides distinct representations for positive and negative
floating-point zeros, (zerop -0.0) is always true. It is an error if the
argument number is not a number.

[Function]
plusp number

This predicate is true if number is strictly greater than zero, and is false
otherwise. It is an error if the argument number is not a non-complex number.

[Function]
minusp number

This predicate is true if number is strictly less than zero, and is false
otherwise. Regardless of whether an implementation provides distinct
representations for positive and negative floating-point zeros, (minusp -0.0)
is always false. (The function float-sign may be used to distinguish a negative
zero.) It is an error if the argument number is not a non-complex number.

[Function]
oddp integer

This predicate is true if the argument integer is odd (not divisible by 2), and
otherwise is false. It is an error if the argument is not an integer.

[Function]
evenp integer

This predicate is true if the argument integer is even (divisible by 2), and
otherwise is false. It is an error if the argument is not an integer.

See also the data-type predicates integerp, rationalp, floatp, complexp, and
numberp.

-------------------------------------------------------------------------------

12.3. Comparisons on Numbers

Each of the functions in this section requires that its arguments all be
numbers; to call one with a non-number is an error. Unless otherwise specified,
each works on all types of numbers, automatically performing any required
coercions when arguments are of different types.

[Function]
= number &rest more-numbers
/= number &rest more-numbers
< number &rest more-numbers
> number &rest more-numbers
<= number &rest more-numbers
>= number &rest more-numbers

These functions each take one or more arguments. If the sequence of arguments
satisfies a certain condition:

=            all the same
/=           all different
<            monotonically increasing
>            monotonically decreasing
<=           monotonically nondecreasing
>=           monotonically nonincreasing

then the predicate is true, and otherwise is false. Complex numbers may be
compared using = and /=, but the others require non-complex arguments. Two
complex numbers are considered equal by = if their real parts are equal and
their imaginary parts are equal according to =. A complex number may be
compared with a non-complex number with = or /=. For example:

(= 3 3) is true.                (/= 3 3) is false.
(= 3 5) is false.               (/= 3 5) is true.
(= 3 3 3 3) is true.            (/= 3 3 3 3) is false.
(= 3 3 5 3) is false.           (/= 3 3 5 3) is false.
(= 3 6 5 2) is false.           (/= 3 6 5 2) is true.
(= 3 2 3) is false.             (/= 3 2 3) is false.
(< 3 5) is true.                (<= 3 5) is true.
(< 3 -5) is false.              (<= 3 -5) is false.
(< 3 3) is false.               (<= 3 3) is true.
(< 0 3 4 6 7) is true.          (<= 0 3 4 6 7) is true.
(< 0 3 4 4 6) is false.         (<= 0 3 4 4 6) is true.
(> 4 3) is true.                (>= 4 3) is true.
(> 4 3 2 1 0) is true.          (>= 4 3 2 1 0) is true.
(> 4 3 3 2 0) is false.         (>= 4 3 3 2 0) is true.
(> 4 3 1 2 0) is false.         (>= 4 3 1 2 0) is false.
(= 3) is true.                  (/= 3) is true.
(< 3) is true.                  (<= 3) is true.
(= 3.0 #C(3.0 0.0)) is true.    (/= 3.0 #C(3.0 1.0)) is true.
(= 3 3.0) is true.              (= 3.0s0 3.0d0) is true.
(= 0.0 -0.0) is true.           (= 5/2 2.5) is true.
(> 0.0 -0.0) is false.          (= 0 -0.0) is true.

With two arguments, these functions perform the usual arithmetic comparison
tests. With three or more arguments, they are useful for range checks, as shown
in the following example:

(<= 0 x 9)                      ;true if x is between 0 and 9, inclusive
(< 0.0 x 1.0)                   ;true if x is between 0.0 and 1.0, exclusive
(< -1 j (length s))             ;true if j is a valid index for s
(<= 0 j k (- (length s) 1))     ;true if j and k are each valid
                                ; indices for s and jk

-------------------------------------------------------------------------------
Rationale: The ``unequality'' relation is called /= rather than <> (the name
used in Pascal) for two reasons. First, /= of more than two arguments is not
the same as the or of < and > of those same arguments. Second, unequality is
meaningful for complex numbers even though < and > are not. For both reasons it
would be misleading to associate unequality with the names of < and >.
-------------------------------------------------------------------------------
Compatibility note: In Common Lisp, the comparison operations perform
``mixed-mode'' comparisons: (= 3 3.0) is true. In MacLisp, there must be
exactly two arguments, and they must be either both fixnums or both
floating-point numbers. To compare two numbers for numerical equality and type
equality, use eql.
-------------------------------------------------------------------------------

[Function]
max number &rest more-numbers
min number &rest more-numbers

The arguments may be any non-complex numbers. max returns the argument that is
greatest (closest to positive infinity). min returns the argument that is least
(closest to negative infinity).

For max, if the arguments are a mixture of rationals and floating-point
numbers, and the largest argument is a rational, then the implementation is
free to produce either that rational or its floating-point approximation; if
the largest argument is a floating-point number of a smaller format than the
largest format of any floating-point argument, then the implementation is free
to return the argument in its given format or expanded to the larger format.
More concisely, the implementation has the choice of returning the largest
argument as is or applying the rules of floating-point contagion, taking all
the arguments into consideration for contagion purposes. Also, if two or more
of the arguments are equal, then any one of them may be chosen as the value to
return. Similar remarks apply to min (replacing ``largest argument'' by
``smallest argument'').

(max 6 12) => 12                (min 6 12) => 6
(max -6 -12) => -6              (min -6 -12) => -12
(max 1 3 2 -7) => 3             (min 1 3 2 -7) => -7
(max -2 3 0 7) => 7             (min -2 3 0 7) => -2
(max 3) => 3                    (min 3) => 3
(max 5.0 2) => 5.0              (min 5.0 2) => 2 or 2.0
(max 3.0 7 1) => 7 or 7.0       (min 3.0 7 1) => 1 or 1.0
(max 1.0s0 7.0d0) => 7.0d0
(min 1.0s0 7.0d0) => 1.0s0 or 1.0d0
(max 3 1 1.0s0 1.0d0) => 3 or 3.0d0
(min 3 1 1.0s0 1.0d0) => 1 or 1.0s0 or 1.0d0

-------------------------------------------------------------------------------

12.4. Arithmetic Operations

Each of the functions in this section requires that its arguments all be
numbers; to call one with a non-number is an error. Unless otherwise specified,
each works on all types of numbers, automatically performing any required
coercions when arguments are of different types.

[Function]
+ &rest numbers

This returns the sum of the arguments. If there are no arguments, the result is
0, which is an identity for this operation.

-------------------------------------------------------------------------------
Compatibility note: While + is compatible with its use in Lisp Machine Lisp, it
is incompatible with MacLisp, which uses + for fixnum-only addition.
-------------------------------------------------------------------------------

[Function]
- number &rest more-numbers

The function -, when given one argument, returns the negative of that argument.

The function -, when given more than one argument, successively subtracts from
the first argument all the others, and returns the result. For example, (- 3 4
5) => -6.

-------------------------------------------------------------------------------
Compatibility note: While - is compatible with its use in Lisp Machine Lisp, it
is incompatible with MacLisp, which uses - for fixnum-only subtraction. Also, -
differs from difference as used in most Lisp systems in the case of one
argument.
-------------------------------------------------------------------------------

[Function]
* &rest numbers

This returns the product of the arguments. If there are no arguments, the
result is 1, which is an identity for this operation.

-------------------------------------------------------------------------------
Compatibility note: While * is compatible with its use in Lisp Machine Lisp, it
is incompatible with MacLisp, which uses * for fixnum-only multiplication.
-------------------------------------------------------------------------------

[Function]
/ number &rest more-numbers

The function /, when given more than one argument, successively divides the
first argument by all the others and returns the result.

[change_begin]
It is generally accepted that it is an error for any argument other than the
first to be zero.
[change_end]

With one argument, / reciprocates the argument.

[change_begin]
It is generally accepted that it is an error in this case for the argument to
be zero.
[change_end]

/ will produce a ratio if the mathematical quotient of two integers is not an
exact integer. For example:

(/ 12 4) => 3
(/ 13 4) => 13/4
(/ -8) => -1/8
(/ 3 4 5) => 3/20

To divide one integer by another producing an integer result, use one of the
functions floor, ceiling, truncate, or round.

If any argument is a floating-point number, then the rules of floating-point
contagion apply.

-------------------------------------------------------------------------------
Compatibility note: What / does is totally unlike what the usual // or quotient
operator does. In most Lisp systems, quotient behaves like / except when
dividing integers, in which case it behaves like truncate of two arguments;
this behavior is mathematically intractable, leading to such anomalies as

(quotient 1.0 2.0) => 0.5   but   (quotient 1 2) => 0

In contrast, the Common Lisp function / produces these results:

(/ 1.0 2.0) => 0.5          and   (/ 1 2) => 1/2

In practice quotient is used only when one is sure that both arguments are
integers, or when one is sure that at least one argument is a floating-point
number. / is tractable for its purpose and works for any numbers.
-------------------------------------------------------------------------------

[Function]
1+ number
1- number

(1+ x) is the same as (+ x 1).

(1- x) is the same as (- x 1). Note that the short name may be confusing: (1-
x) does not mean 1-x; rather, it means x-1.

-------------------------------------------------------------------------------
Rationale: These are included primarily for compatibility with MacLisp and Lisp
Machine Lisp. Some programmers prefer always to write (+ x 1) and (- x 1)
instead of (1+ x) and (1- x).
-------------------------------------------------------------------------------
Implementation note: Compiler writers are very strongly encouraged to ensure
that (1+ x) and (+ x 1) compile into identical code, and similarly for (1- x)
and (- x 1), to avoid pressure on a Lisp programmer to write possibly less
clear code for the sake of efficiency. This can easily be done as a
source-language transformation.
-------------------------------------------------------------------------------

[Macro]
incf place [delta]
decf place [delta]

The number produced by the form delta is added to (incf) or subtracted from
(decf) the number in the generalized variable named by place, and the sum is
stored back into place and returned. The form place may be any form acceptable
as a generalized variable to setf. If delta is not supplied, then the number in
place is changed by 1. For example:

(setq n 0)
(incf n) => 1                   and now n => 1
(decf n 3) => -2                and now n => -2
(decf n -5) => 3                and now n => 3
(decf n) => 2                   and now n => 2

The effect of (incf place delta) is roughly equivalent to

(setf place (+ place delta))

except that the latter would evaluate any subforms of place twice, whereas incf
takes care to evaluate them only once. Moreover, for certain place forms incf
may be significantly more efficient than the setf version.

[change_begin]
X3J13 voted in March 1988 (PUSH-EVALUATION-ORDER)   to clarify order of
evaluation (see section 7.2).
[change_end]

[Function]
conjugate number

This returns the complex conjugate of number. The conjugate of a non-complex
number is itself. For a complex number z,

(conjugate z) == (complex (realpart z) (- (imagpart z)))

For example:

(conjugate #C(3/5 4/5)) => #C(3/5 -4/5)
(conjugate #C(0.0D0 -1.0D0)) => #C(0.0D0 1.0D0)
(conjugate 3.7) => 3.7

[Function]
gcd &rest integers

This returns the greatest common divisor of all the arguments, which must be
integers. The result of gcd is always a non-negative integer. If one argument
is given, its absolute value is returned. If no arguments are given, gcd
returns 0, which is an identity for this operation. For three or more
arguments,

(gcd a b c ... z) == (gcd (gcd a b) c ... z)

Here are some examples of the use of gcd:

(gcd 91 -49) => 7
(gcd 63 -42 35) => 7
(gcd 5) => 5
(gcd -4) => 4
(gcd) => 0

[Function]
lcm integer &rest more-integers

This returns the least common multiple of its arguments, which must be
integers. The result of lcm is always a non-negative integer. For two arguments
that are not both zero,

(lcm a b) == (/ (abs (* a b)) (gcd a b))

If one or both arguments are zero,

(lcm a 0) == (lcm 0 a) == 0

For one argument, lcm returns the absolute value of that argument. For three or
more arguments,

(lcm a b c ... z) == (lcm (lcm a b) c ... z)

Some examples:

(lcm 14 35) => 70
(lcm 0 5) => 0
(lcm 1 2 3 4 5 6) => 60

Mathematically, (lcm) should return infinity. Because Common Lisp does not have
a representation for infinity, lcm, unlike gcd, always requires at least one
argument.

[change_begin]
X3J13 voted in January 1989 (LCM-NO-ARGUMENTS)   to specify that (lcm) => 1.

This is one of my biggest boners. The identity for lcm is of course 1, not
infinity, and so (lcm) ought to have been defined to return 1. Sorry about
that, though in point of fact very few users have complained to me that this
mistake in the first edition has cramped their programming style.
[change_end]

-------------------------------------------------------------------------------

12.5. Irrational and Transcendental Functions

Common Lisp provides no data type that can accurately represent irrational
numerical values. The functions in this section are described as if the results
were mathematically accurate, but actually they all produce floating-point
approximations to the true mathematical result in the general case. In some
places mathematical identities are set forth that are intended to elucidate the
meanings of the functions; however, two mathematically identical expressions
may be computationally different because of errors inherent in the
floating-point approximation process.

When the arguments to a function in this section are all rational and the true
mathematical result is also (mathematically) rational, then unless otherwise
noted an implementation is free to return either an accurate result of type
rational or a single-precision floating-point approximation. If the arguments
are all rational but the result cannot be expressed as a rational number, then
a single-precision floating-point approximation is always returned.

[change_begin]
X3J13 voted in March 1989 (COMPLEX-RATIONAL-RESULT)   to clarify that the
provisions of the previous paragraph apply to complex numbers. If the arguments
to a function are all of type (or rational (complex rational)) and the true
mathematical result is (mathematically) a complex number with rational real and
imaginary parts, then unless otherwise noted an implementation is free to
return either an accurate result of type (or rational (complex rational)) or a
single-precision floating-point approximation of type single-float (permissible
only if the imaginary part of the true mathematical result is zero) or (complex
single-float). If the arguments are all of type (or rational (complex
rational)) but the result cannot be expressed as a rational or complex rational
number, then the returned value will be of type single-float (permissible only
if the imaginary part of the true mathematical result is zero) or (complex
single-float).
[change_end]

The rules of floating-point contagion and complex contagion are effectively
obeyed by all the functions in this section except expt, which treats some
cases of rational exponents specially. When, possibly after contagious
conversion, all of the arguments are of the same floating-point or complex
floating-point type, then the result will be of that same type unless otherwise
noted.

-------------------------------------------------------------------------------
Implementation note: There is a ``floating-point cookbook'' by Cody and Waite
[14] that may be a useful aid in implementing the functions defined in this
section.
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

   *  Exponential and Logarithmic Functions
   *  Trigonometric and Related Functions
   *  Branch Cuts, Principal Values, and Boundary Conditions in the Complex
     Plane

-------------------------------------------------------------------------------

12.5.1. Exponential and Logarithmic Functions

Along with the usual one-argument and two-argument exponential and logarithm
functions, sqrt is considered to be an exponential function, because it raises
a number to the power 1/2.

[Function]
exp number

Returns e raised to the power number, where e is the base of the natural
logarithms.

[Function]
expt base-number power-number

Returns base-number raised to the power power-number. If the base-number is of
type rational and the power-number is an integer, the calculation will be exact
and the result will be of type rational; otherwise a floating-point
approximation may result.

[change_begin]
X3J13 voted in March 1989 (COMPLEX-RATIONAL-RESULT)   to clarify that
provisions similar to those of the previous paragraph apply to complex numbers.
If the base-number is of type (complex rational) and the power-number is an
integer, the calculation will also be exact and the result will be of type (or
rational (complex rational)); otherwise a floating-point or complex
floating-point approximation may result.
[change_end]

When power-number is 0 (a zero of type integer), then the result is always the
value 1 in the type of base-number, even if the base-number is zero (of any
type). That is:

(expt x 0) == (coerce 1 (type-of x))

If the power-number is a zero of any other data type, then the result is also
the value 1, in the type of the arguments after the application of the
contagion rules, with one exception: it is an error if base-number is zero when
the power-number is a zero not of type integer.

Implementations of expt are permitted to use different algorithms for the cases
of a rational power-number and a floating-point power-number; the motivation is
that in many cases greater accuracy can be achieved for the case of a rational
power-number. For example, (expt pi 16) and (expt pi 16.0) may yield slightly
different results if the first case is computed by repeated squaring and the
second by the use of logarithms. Similarly, an implementation might choose to
compute (expt x 3/2) as if it had been written (sqrt (expt x 3)), perhaps
producing a more accurate result than would (expt x 1.5). It is left to the
implementor to determine the best strategies.

[change_begin]
X3J13 voted in January 1989 (EXPT-RATIO)   to clarify that the preceding remark
is in error, because (sqrt (expt x 3)) does not produce the same value as (expt
x 3/2) in most cases, and to specify that the specification of the principal
value of expt as given in section 12.5.3 should be regarded as definitive.

As an example of the difficulty, let   . Then   , but   . Another example is
x=-1; then   , but   .
[change_end]

The result of expt can be a complex number, even when neither argument is
complex, if base-number is negative and power-number is not an integer. The
result is always the principal complex value. Note that (expt -8 1/3) is not
permitted to return -2; while -2 is indeed one of the cube roots of -8, it is
not the principal cube root, which is a complex number approximately equal to
#C(1.0 1.73205).

[change_begin]
Notice of correction. The first edition gave the incorrect value #C(0.5
1.73205) for the principal cube root of -8. The correct value is #C(1.0
1.73205), that is, 1+SQRT(3)i. I simply don't know what I was thinking of!
[change_end]

[Function]
log number &optional base

Returns the logarithm of number in the base base, which defaults to e, the base
of the natural logarithms. For example:

(log 8.0 2) => 3.0
(log 100.0 10) => 2.0

The result of (log 8 2) may be either 3 or 3.0, depending on the
implementation.

Note that log may return a complex result when given a non-complex argument if
the argument is negative. For example:

(log -1.0) == (complex 0.0 (float pi 0.0))

[change_begin]
X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
floating-point behavior when minus zero is supported. As a part of that vote it
approved a mathematical definition of complex logarithm in terms of real
logarithm, absolute value, arc tangent of two real arguments, and the phase
function as


Logarithm               log|z| + i phase z

This specifies the branch cuts precisely whether minus zero is supported or
not; see phase and atan.
[change_end]

[Function]
sqrt number

Returns the principal square root of number. If the number is not complex but
is negative, then the result will be a complex number. For example:

(sqrt 9.0) => 3.0
(sqrt -9.0) => #c(0.0 3.0)

The result of (sqrt 9) may be either 3 or 3.0, depending on the implementation.
The result of (sqrt -9) may be either #c(0 3) or #c(0.0 3.0).

[change_begin]
X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
floating-point behavior when minus zero is supported. As a part of that vote it
approved a mathematical definition of complex square root in terms of complex
logarithm and exponential functions as


Square root

This specifies the branch cuts precisely whether minus zero is supported or
not; see phase and atan.
[change_end]

[Function]
isqrt integer

Integer square root: the argument must be a non-negative integer, and the
result is the greatest integer less than or equal to the exact positive square
root of the argument. For example:

(isqrt 9) => 3
(isqrt 12) => 3
(isqrt 300) => 17
(isqrt 325) => 18

-------------------------------------------------------------------------------

12.5.2. Trigonometric and Related Functions

Some of the functions in this section, such as abs and signum, are apparently
unrelated to trigonometric functions when considered as functions of real
numbers only. The way in which they are extended to operate on complex numbers
makes the trigonometric connection clear.

[Function]
abs number

Returns the absolute value of the argument. For a non-complex number x,

(abs x) == (if (minusp x) (- x) x)

and the result is always of the same type as the argument.

For a complex number z, the absolute value may be computed as

(sqrt (+ (expt (realpart z) 2) (expt (imagpart z) 2)))

-------------------------------------------------------------------------------
Implementation note: The careful implementor will not use this formula directly
for all complex numbers but will instead handle very large or very small
components specially to avoid intermediate overflow or underflow.
-------------------------------------------------------------------------------

For example:

(abs #c(3.0 -4.0)) => 5.0

The result of (abs #c(3 4)) may be either 5 or 5.0, depending on the
implementation.

[Function]
phase number

The phase of a number is the angle part of its polar representation as a
complex number. That is,

(phase z) == (atan (imagpart z) (realpart z))

[old_change_begin]
The result is in radians, in the range -pi (exclusive) to +pi (inclusive). The
phase of a positive non-complex number is zero; that of a negative non-complex
number is pi. The phase of zero is arbitrarily defined to be zero.
[old_change_end]

[change_begin]
X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
floating-point behavior when minus zero is supported; phase is still defined in
terms of atan as above, but thanks to a change in atan the range of phase
becomes -pi inclusive to +pi inclusive. The value - results from an argument
whose real part is negative and whose imaginary part is minus zero. The phase
function therefore has a branch cut along the negative real axis. The phase of
+0+0i is +0, of +0-0i is -0, of -0+0i is +pi, and of -0-0i is -pi.
[change_end]

If the argument is a complex floating-point number, the result is a
floating-point number of the same type as the components of the argument. If
the argument is a floating-point number, the result is a floating-point number
of the same type. If the argument is a rational number or complex rational
number, the result is a single-format floating-point number.

[Function]
signum number

By definition,

(signum x) == (if (zerop x) x (/ x (abs x)))

For a rational number, signum will return one of -1, 0, or 1 according to
whether the number is negative, zero, or positive. For a floating-point number,
the result will be a floating-point number of the same format whose value is
-1, 0, or 1. For a complex number z, (signum z) is a complex number of the same
phase but with unit magnitude, unless z is a complex zero, in which case the
result is z. For example:

(signum 0) => 0
(signum -3.7L5) => -1.0L0
(signum 4/5) => 1
(signum #C(7.5 10.0)) => #C(0.6 0.8)
(signum #C(0.0 -14.7)) => #C(0.0 -1.0)

For non-complex rational numbers, signum is a rational function, but it may be
irrational for complex arguments.

[Function]
sin radians
cos radians
tan radians

sin returns the sine of the argument, cos the cosine, and tan the tangent. The
argument is in radians. The argument may be complex.

[Function]
cis radians

This computes   . The name cis means ``cos + i sin,'' because   . The argument
is in radians and may be any non-complex number. The result is a complex number
whose real part is the cosine of the argument and whose imaginary part is the
sine. Put another way, the result is a complex number whose phase is the equal
to the argument (mod 2) and whose magnitude is unity.

-------------------------------------------------------------------------------
Implementation note: Often it is cheaper to calculate the sine and cosine of a
single angle together than to perform two disjoint calculations.
-------------------------------------------------------------------------------

[Function]
asin number
acos number

asin returns the arc sine of the argument, and acos the arc cosine. The result
is in radians. The argument may be complex.

The arc sine and arc cosine functions may be defined mathematically for an
argument z as follows:


Arc sine

Arc cosine

Note that the result of asin or acos may be complex even if the argument is not
complex; this occurs when the absolute value of the argument is greater than 1.

[change_begin]
Kahan [25] suggests for acos the defining formula


Arc cosine

or even the much simpler   . Both equations are mathematically equivalent to
the formula shown above.
[change_end]

-------------------------------------------------------------------------------
Implementation note: These formulae are mathematically correct, assuming
completely accurate computation. They may be terrible methods for
floating-point computation. Implementors should consult a good text on
numerical analysis. The formulae given above are not necessarily the simplest
ones for real-valued computations, either; they are chosen to define the branch
cuts in desirable ways for the complex case.
-------------------------------------------------------------------------------

[Function]
atan y &optional x

An arc tangent is calculated and the result is returned in radians.

With two arguments y and x, neither argument may be complex. The result is the
arc tangent of the quantity y/x. The signs of y and x are used to derive
quadrant information; moreover, x may be zero provided y is not zero. The value
of atan is always between -pi (exclusive) and +pi (inclusive). The following
table details various special cases.



[change_begin]
X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
floating-point behavior when minus zero is supported. When there is a minus
zero, the preceding table must be modified slightly:



Note that the case y=0,x=0 is an error in the absence of minus zero, but the
four cases y=0,x=0 are defined in the presence of minus zero.
[change_end]

[old_change_begin]
With only one argument y, the argument may be complex. The result is the arc
tangent of y, which may be defined by the following formula:

Arc tangent

[old_change_end]

-------------------------------------------------------------------------------
Implementation note: This formula is mathematically correct, assuming
completely accurate computation. It may be a terrible method for floating-point
computation. Implementors should consult a good text on numerical analysis. The
formula given above is not necessarily the simplest one for real-valued
computations, either; it is chosen to define the branch cuts in desirable ways
for the complex case.
-------------------------------------------------------------------------------

[change_begin]
X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to replace the
preceding formula with the formula

                        log(1+iy) - log(1-iy)
Arc tangent             ---------------------
                                  2i

This change alters the direction of continuity for the branch cuts, which
alters the result returned by atan only for arguments on the imaginary axis
that are of magnitude greater than 1. See section 12.5.3 for further details.
[change_end]

For a non-complex argument y, the result is non-complex and lies between    and
   (both exclusive).

-------------------------------------------------------------------------------
Compatibility note: MacLisp has a function called atan whose range is from 0 to
2. Almost every other programming language (ANSI Fortran, IBM PL/1, Interlisp)
has a two-argument arc tangent function with range -pi to +pi. Lisp Machine
Lisp provides two two-argument arc tangent functions, atan (compatible with
MacLisp) and atan2 (compatible with all others).

Common Lisp makes two-argument atan the standard one with range -pi to +pi.
Observe that this makes the one-argument and two-argument versions of atan
compatible in the sense that the branch cuts do not fall in different places.
The Interlisp one-argument function arctan has a range from 0 to pi, while
nearly every other programming language provides the range    to    for
one-argument arc tangent! Nevertheless, since Interlisp uses the standard
two-argument version of arc tangent, its branch cuts are inconsistent anyway.
-------------------------------------------------------------------------------

[Constant]
pi

This global variable has as its value the best possible approximation to pi in
long floating-point format. For example:

(defun sind (x)     ;The argument is in degrees
  (sin (* x (/ (float pi x) 180))))

An approximation to pi in some other precision can be obtained by writing
(float pi x), where x is a floating-point number of the desired precision, or
by writing (coerce pi type), where type is the name of the desired type, such
as short-float.

[Function]
sinh number
cosh number
tanh number
asinh number
acosh number
atanh number

[old_change_begin]
These functions compute the hyperbolic sine, cosine, tangent, arc sine, arc
cosine, and arc tangent functions, which are mathematically defined for an
argument z as follows:

Hyperbolic sine
Hyperbolic cosine
Hyperbolic tangent
Hyperbolic arc sine
Hyperbolic arc cosine
Hyperbolic arc tangent           WRONG!

[old_change_end]

[change_begin]
WARNING! The formula shown above for hyperbolic arc tangent is incorrect. It is
not a matter of incorrect branch cuts; it simply does not compute anything like
a hyperbolic arc tangent. This unfortunate error in the first edition was the
result of mistranscribing a (correct) APL formula from Penfield's paper [36].
The formula should have been transcribed as


Hyperbolic arc tangent

A proposal was submitted to X3J13 in September 1989 to replace the formulae for
acosh and atanh. See section 12.5.3 for further discussion.
[change_end]

Note that the result of acosh may be complex even if the argument is not
complex; this occurs when the argument is less than 1. Also, the result of
atanh may be complex even if the argument is not complex; this occurs when the
absolute value of the argument is greater than 1.

-------------------------------------------------------------------------------
Implementation note: These formulae are mathematically correct, assuming
completely accurate computation. They may be terrible methods for
floating-point computation. Implementors should consult a good text on
numerical analysis. The formulae given above are not necessarily the simplest
ones for real-valued computations, either; they are chosen to define the branch
cuts in desirable ways for the complex case.
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

12.5.3. Branch Cuts, Principal Values, and Boundary Conditions in the Complex
Plane

Many of the irrational and transcendental functions are multiply defined in the
complex domain; for example, there are in general an infinite number of complex
values for the logarithm function. In each such case, a principal value must be
chosen for the function to return. In general, such values cannot be chosen so
as to make the range continuous; lines in the domain called branch cuts must be
defined, which in turn define the discontinuities in the range.

Common Lisp defines the branch cuts, principal values, and boundary conditions
for the complex functions following a proposal for complex functions in APL
[36]. The contents of this section are borrowed largely from that proposal.

-------------------------------------------------------------------------------
Compatibility note: The branch cuts defined here differ in a few very minor
respects from those advanced by W. Kahan, who considers not only the ``usual''
definitions but also the special modifications necessary for IEEE proposed
floating-point arithmetic, which has infinities and minus zero as explicit
computational objects. For example, he proposes that SQRT(-4+0i)=2i, but
SQRT(-4-0i)=-2i.

It may be that the differences between the APL proposal and Kahan's proposal
will be ironed out. If so, Common Lisp may be changed as necessary to be
compatible with these other groups. Any changes from the specification below
are likely to be quite minor, probably concerning primarily questions of which
side of a branch cut is continuous with the cut itself.
-------------------------------------------------------------------------------

[change_begin]
Indeed, X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to alter the
direction of continuity for the branch cuts of atan, and also
(IEEE-ATAN-BRANCH-CUT)   to address the treatment of branch cuts in
implementations that have a distinct floating-point minus zero.

The treatment of minus zero centers in two-argument atan. If there is no minus
zero, then the branch cut runs just below the negative real axis as before, and
the range of two-argument atan is   . If there is a minus zero, however, then
the branch cut runs precisely on the negative real axis, skittering between
pairs of numbers of the form -x+/-0i, and the range of two-argument atan is   .

The treatment of minus zero by all other irrational and transcendental
functions is then specified by defining those functions in terms of
two-argument atan. First, phase is defined in terms of two-argument atan, and
complex abs in terms of real sqrt; then complex log is defined in terms of
phase, abs, and real log; then complex sqrt in terms of complex log; and
finally all others are defined in terms of these.

Kahan [25] treats these matters in some detail and also suggests specific
algorithms for implementing irrational and transcendental functions in IEEE
standard floating-point arithmetic [23].

Remarks in the first edition about the direction of the continuity of branch
cuts continue to hold in the absence of minus zero and may be ignored if minus
zero is supported; since all branch cuts happen to run along the principal
axes, they run between plus zero and minus zero, and so each sort of zero is
associated with the obvious quadrant.
[change_end]

sqrt
     The branch cut for square root lies along the negative real axis,
     continuous with quadrant II. The range consists of the right half-plane,
     including the non-negative imaginary axis and excluding the negative
     imaginary axis.

[change_begin]

     X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
     floating-point behavior when minus zero is supported. As a part of that
     vote it approved a mathematical definition of complex square root:



     This defines the branch cuts precisely, whether minus zero is supported or
     not.

[change_end]

phase
     The branch cut for the phase function lies along the negative real axis,
     continuous with quadrant II. The range consists of that portion of the
     real axis between -pi (exclusive) and pi (inclusive).

[change_begin]

     X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
     floating-point behavior when minus zero is supported. As a part of that
     vote it approved a mathematical definition of phase:



     where Fz is the imaginary part of z and Rz the real part of z. This
     defines the branch cuts precisely, whether minus zero is supported or not.

[change_end]

log  The branch cut for the logarithm function of one argument (natural
     logarithm) lies along the negative real axis, continuous with quadrant II.
     The domain excludes the origin. For a complex number z, z is defined to be

      log z = (log|z|) + i(phase z)

     Therefore the range of the one-argument logarithm function is that strip
     of the complex plane containing numbers with imaginary parts between -pi
     (exclusive) and pi (inclusive).

     [change_begin]
     The X3J13 vote on minus zero (IEEE-ATAN-BRANCH-CUT)   would alter that
     exclusive bound of -pi to be inclusive if minus zero is supported.
     [change_end]

     The two-argument logarithm function is defined as   . This defines the
     principal values precisely. The range of the two-argument logarithm
     function is the entire complex plane. It is an error if z is zero. If z is
     non-zero and b is zero, the logarithm is taken to be zero.

exp  The simple exponential function has no branch cut.

expt
     The two-argument exponential function is defined as   . This defines the
     principal values precisely. The range of the two-argument exponential
     function is the entire complex plane. Regarded as a function of x, with b
     fixed, there is no branch cut. Regarded as a function of b, with x fixed,
     there is in general a branch cut along the negative real axis, continuous
     with quadrant II. The domain excludes the origin. By definition,   . If
     b=0 and the real part of x is strictly positive, then   . For all other
     values of x,    is an error.

asin
     The following definition for arc sine determines the range and branch
     cuts:



     This is equivalent to the formula



     recommended by Kahan [25].

     The branch cut for the arc sine function is in two pieces: one along the
     negative real axis to the left of -1 (inclusive), continuous with quadrant
     II, and one along the positive real axis to the right of 1 (inclusive),
     continuous with quadrant IV. The range is that strip of the complex plane
     containing numbers whose real part is between    and   . A number with
     real part equal to    is in the range if and only if its imaginary part is
     non-negative; a number with real part equal to    is in the range if and
     only if its imaginary part is non-positive.

acos
     The following definition for arc cosine determines the range and branch
     cuts:



     or, which is equivalent,



     The branch cut for the arc cosine function is in two pieces: one along the
     negative real axis to the left of -1 (inclusive), continuous with quadrant
     II, and one along the positive real axis to the right of 1 (inclusive),
     continuous with quadrant IV. This is the same branch cut as for arc sine.
     The range is that strip of the complex plane containing numbers whose real
     part is between zero and pi. A number with real part equal to zero is in
     the range if and only if its imaginary part is non-negative; a number with
     real part equal to pi is in the range if and only if its imaginary part is
     non-positive.

atan
     The following definition for (one-argument) arc tangent determines the
     range and branch cuts:



[old_change_begin]

     Beware of simplifying this formula; ``obvious'' simplifications are likely
     to alter the branch cuts or the values on the branch cuts incorrectly.

     The branch cut for the arc tangent function is in two pieces: one along
     the positive imaginary axis above i (exclusive), continuous with quadrant
     II, and one along the negative imaginary axis below -i (exclusive),
     continuous with quadrant IV. The points i and -i are excluded from the
     domain. The range is that strip of the complex plane containing numbers
     whose real part is between    and   . A number with real part equal to
     is in the range if and only if its imaginary part is strictly positive; a
     number with real part equal to    is in the range if and only if its
     imaginary part is strictly negative. Thus the range of the arc tangent
     function is identical to that of the arc sine function with the points
     and    excluded.

[old_change_end]

[change_begin]

     X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to replace the
     formula shown above with the formula



     This is equivalent to the formula



     recommended by Kahan [25]. It causes the upper branch cut to be continuous
     with quadrant I rather than quadrant II, and the lower branch cut to be
     continuous with quadrant III rather than quadrant IV; otherwise it agrees
     with the formula of the first edition. Therefore this change alters the
     result returned by atan only for arguments on the positive imaginary axis
     that are of magnitude greater than 1. The full description for this new
     formula is as follows.

     The branch cut for the arc tangent function is in two pieces: one along
     the positive imaginary axis above i (exclusive), continuous with quadrant
     I, and one along the negative imaginary axis below -i (exclusive),
     continuous with quadrant III. The points i and -i are excluded from the
     domain. The range is that strip of the complex plane containing numbers
     whose real part is between    and   . A number with real part equal to
     is in the range if and only if its imaginary part is strictly negative; a
     number with real part equal to    is in the range if and only if its
     imaginary part is strictly positive. Thus the range of the arc tangent
     function is not identical to that of the arc sine function.

[change_end]

asinh
     The following definition for the inverse hyperbolic sine determines the
     range and branch cuts:



     The branch cut for the inverse hyperbolic sine function is in two pieces:
     one along the positive imaginary axis above i (inclusive), continuous with
     quadrant I, and one along the negative imaginary axis below -i
     (inclusive), continuous with quadrant III. The range is that strip of the
     complex plane containing numbers whose imaginary part is between    and
       . A number with imaginary part equal to    is in the range if and only
     if its real part is non-positive; a number with imaginary part equal to
     is in the range if and only if its real part is non-negative.

acosh
     The following definition for the inverse hyperbolic cosine determines the
     range and branch cuts:



[change_begin]

     Kahan [25] suggests the formula



     pointing out that it yields the same principal value but eliminates a
     gratuitous removable singularity at z=-1. A proposal was submitted to
     X3J13 in September 1989 to replace the formula acosh with that recommended
     by Kahan. There is a good possibility that it will be adopted.

[change_end]

     The branch cut for the inverse hyperbolic cosine function lies along the
     real axis to the left of 1 (inclusive), extending indefinitely along the
     negative real axis, continuous with quadrant II and (between 0 and 1) with
     quadrant I. The range is that half-strip of the complex plane containing
     numbers whose real part is non-negative and whose imaginary part is
     between -pi (exclusive) and pi (inclusive). A number with real part zero
     is in the range if its imaginary part is between zero (inclusive) and pi
     (inclusive).

atanh
     The following definition for the inverse hyperbolic tangent determines the
     range and branch cuts:

[old_change_begin]

             WRONG!

[old_change_end]

[change_begin]

     WARNING! The formula shown above for hyperbolic arc tangent is incorrect.
     It is not a matter of incorrect branch cuts; it simply does not compute
     anything like a hyperbolic arc tangent. This unfortunate error in the
     first edition was the result of mistranscribing a (correct) APL formula
     from Penfield's paper [36]. The formula should have been transcribed as



[change_end]

[old_change_begin]

     Beware of simplifying this formula; ``obvious'' simplifications are likely
     to alter the branch cuts or the values on the branch cuts incorrectly.

     The branch cut for the inverse hyperbolic tangent function is in two
     pieces: one along the negative real axis to the left of -1 (inclusive),
     continuous with quadrant III, and one along the positive real axis to the
     right of 1 (inclusive), continuous with quadrant I. The points -1 and 1
     are excluded from the domain. The range is that strip of the complex plane
     containing numbers whose imaginary part is between    and   . A number
     with imaginary part equal to    is in the range if and only if its real
     part is strictly negative; a number with imaginary part equal to    is in
     the range if and only if its real part is strictly positive. Thus the
     range of the inverse hyperbolic tangent function is identical to that of
     the inverse hyperbolic sine function with the points    and    excluded.

[old_change_end]

[change_begin]

     A proposal was submitted to X3J13 in September 1989 to replace the formula
     atanh with that recommended by Kahan [25]:



     There is a good possibility that it will be adopted. If it is, the
     complete description of the branch cuts of atanh will then be as follows.

     The branch cut for the inverse hyperbolic tangent function is in two
     pieces: one along the negative real axis to the left of -1 (inclusive),
     continuous with quadrant II, and one along the positive real axis to the
     right of 1 (inclusive), continuous with quadrant IV. The points -1 and 1
     are excluded from the domain. The range is that strip of the complex plane
     containing numbers whose imaginary part is between    and   . A number
     with imaginary part equal to    is in the range if and only if its real
     part is strictly positive; a number with imaginary part equal to    is in
     the range if and only if its real part is strictly negative. Thus the
     range of the inverse hyperbolic tangent function is not the same as that
     of the inverse hyperbolic sine function.

[change_end]

With these definitions, the following useful identities are obeyed throughout
the applicable portion of the complex domain, even on the branch cuts:

sin iz = i sinh z  sinh iz = i sin z        arctan iz = i arctanh z
cos iz = cosh z    cosh iz = cos z          arcsinh iz = i arcsin z
tan iz = i tanh z  arcsin iz = i arcsinh z  arctanh iz = i arctan z

[change_begin]
I thought it would be useful to provide some graphs illustrating the behavior
of the irrational and transcendental functions in the complex plane. It also
provides an opportunity to show off the Common Lisp code that was used to
generate them.

Imagine the complex plane to be decorated as follows. The real and imaginary
axes are painted with thick lines. Parallels from the axes on both sides at
distances of 1, 2, and 3 are painted with thin lines; these parallels are
doubly infinite lines, as are the axes. Four annuli (rings) are painted in
gradated shades of gray. Ring 1, the inner ring, consists of points whose
radial distances from the origin lie in the range [1/4, 1/2]; ring 2 is in the
radial range [3/4, 1]; ring 3, in the range [pi/2, 2]; and ring 4, in the range
[3, pi]. Ring j is divided into    equal sectors, with each sector painted a
different shade of gray, darkening as one proceeds counterclockwise from the
positive real axis.

We can illustrate the behavior of a numerical function f by considering how it
maps the complex plane to itself. More specifically, consider each point z of
the decorated plane. We decorate a new plane by coloring the point f(z) with
the same color that point z had in the original decorated plane. In other
words, the newly decorated plane illustrates how the f maps the axes, other
horizontal and vertical lines, and annuli.

In each figure we will show only a fragment of the complex plane, with the real
axis horizontal in the usual manner (-pi to the left, +pi to the right) and the
imaginary axis vertical (-i below, +pii above). Each fragment shows a region
containing points whose real and imaginary parts are in the range [-4.1, 4.1].
The axes of the new plane are shown as very thin lines, with large tick marks
at integer coordinates and somewhat smaller tick marks at multiples of   .

Figure 12-1 shows the result of plotting the identity function (quite
literally); the graph exhibits the decoration of the original plane.

Figures 12-2 through 12-20 show the graphs for the functions sqrt, exp, log,
sin, asin, cos, acos, tan, atan, sinh, asinh, cosh, acosh, tanh, and atanh, and
as a bonus, the graphs for the functions   ,   ,   , and   . All of these are
related to the trigonometric functions in various ways. For example, if   ,
then   , and if   , then   . It is instructive to examine the graph for    and
try to visualize how it transforms the graph for sin into the graph for cos.

Each figure is accompanied by a commentary on what maps to what and other
interesting features. None of this material is terribly new; much of it may be
found in any good textbook on complex analysis. I believe that the particular
form in which the graphs are presented is novel, as well as the fact that the
graphs have been generated as PostScript [1] code by Common Lisp code. This
PostScript code was then fed directly to the typesetting equipment that set the
pages for this book. Samples of this PostScript code follow the figures
themselves, after which the code for the entire program is presented.

In the commentaries that accompany the figures I sometimes speak of mapping the
points +/- infinity or +/- infinity i. When I say that function f maps
+infinity to a certain point z, I mean that



Similarly, when I say that f maps -i to z, I mean that



In other words, I am considering a limit as one travels out along one of the
main axes. I also speak in a similar manner of mapping to one of these
infinities.

-------------------------------------------------------------------------------

   *   Identity Plot
   *   Square Root Function
   *  Exponential Function
   *  Natural Logarithm Function
   *  Function (z-1)/(z+1)
   *  Function (1+z)/(1-z)
   *  Sine Function
   *  Arc Sine Function
   *  Cosine Function
   *  Arc Cosine Function
   *  Tanget Function
   *  Arc Tangent Function
   *  Hyperbolic Sine Function
   *  Hyperbolic Arc Sine Function
   *  Hyperbolic Cosine Function
   *  Hyperbolic Arc Cosine Function
   *  Hyperbolic Tangent Function
   *   Hyperbolic Arc Tangent Function
   *  SQRT(1 - SQR(z))
   *   SQRT(1 + SQR(z))

Here is a sample of the PostScript code that generated figure 12-1, showing the
initial scaling, translation, and clipping parameters; the code for one sector
of the innermost annulus; and the code for the negative imaginary axis. Comment
lines indicate how path or boundary segments were generated separately and then
spliced (in order to allow for the places that a singularity might lurk, in
which case the generating code can ``inch up'' to the problematical argument
value).

The size of the entire PostScript file for the identity function was about 68
kilobytes (2757 lines, including comments). The smallest files were the plots
for atan and atanh, about 65 kilobytes apiece; the largest were the plots for
sin, cos, sinh, and cosh, about 138 kilobytes apiece.

% PostScript file for plot of function IDENTITY
% Plot is to fit in a region 4.666666666666667 inches square
%  showing axes extending 4.1 units from the origin.

40.97560975609756 40.97560975609756 scale
4.1 4.1 translate
newpath
  -4.1 -4.1 moveto
  4.1 -4.1 lineto
  4.1 4.1 lineto
  -4.1 4.1 lineto
  closepath
clip
% Moby grid for function IDENTITY
% Annulus 0.25 0.5 4 0.97 0.45
% Sector from 4.7124 to 6.2832 (quadrant 3)
newpath
  0.0 -0.25 moveto
  0.0 -0.375 lineto
  %middle radial
  0.0 -0.375 lineto
  0.0 -0.5 lineto
  %end radial
  0.0 -0.5 lineto
  0.092 -0.4915 lineto
  0.1843 -0.4648 lineto
  0.273 -0.4189 lineto
  0.3536 -0.3536 lineto
  %middle circumferential
  0.3536 -0.3536 lineto
  0.413 -0.2818 lineto
  0.4594 -0.1974 lineto
  0.4894 -0.1024 lineto
  0.5 0.0 lineto
  %end circumferential
  0.5 0.0 lineto
  0.375 0.0 lineto
  %middle radial
  0.375 0.0 lineto
  0.25 0.0 lineto
  %end radial
  0.25 0.0 lineto
  0.2297 -0.0987 lineto
  0.1768 -0.1768 lineto
  %middle circumferential
  0.1768 -0.1768 lineto
  0.0922 -0.2324 lineto
  0.0 -0.25 lineto
  %end circumferential
  closepath
currentgray   0.45 setgray   fill   setgray

[2598 lines omitted]

% Vertical line from (0.0, -0.5) to (0.0, 0.0)
newpath
  0.0 -0.5 moveto
  0.0 0.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -0.5) to (0.0, -1.0)
newpath
  0.0 -0.5 moveto
  0.0 -1.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -2.0) to (0.0, -1.0)
newpath
  0.0 -2.0 moveto
  0.0 -1.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -2.0) to (0.0, -1.1579208923731617E77)
newpath
  0.0 -2.0 moveto
  0.0 -6.3553 lineto
  0.0 -6.378103166302659 lineto
  0.0 -6.378103166302659 lineto
  0.0 -6.378103166302659 lineto
0.05 setlinewidth   1 setlinecap  stroke

[84 lines omitted]

% End of PostScript file for plot of function IDENTITY

Here is the program that generated the PostScript code for the graphs shown in
figures 12-1 through 12-20. It contains a mixture of fairly general mechanisms
and ad hoc kludges for plotting functions of a single complex argument while
gracefully handling extremely large and small values, branch cuts,
singularities, and periodic behavior. The aim was to provide a simple user
interface that would not require the caller to provide special advice for each
function to be plotted. The file for figure 12-1, for example, was generated by
the call (picture 'identity), which resulted in the writing of a file named
identity-plot.ps.

The program assumes that any periodic behavior will have a period that is a
multiple of 2pi; that branch cuts will fall along the real or imaginary axis;
and that singularities or very large or small values will occur only at the
origin, at 1 or +/-i, or on the boundaries of the annuli (particularly those
with radius    or pi). The central function is parametric-path, which accepts
four arguments: two real numbers that are the endpoints of an interval of real
numbers, a function that maps this interval into a path in the complex plane,
and the function to be plotted; the task of parametric-path is to generate
PostScript code (a series of lineto operations) that will plot an approximation
to the image of the parametric path as transformed by the function to be
plotted. Each of the functions hline, vline, -hline, -vline, radial, and
circumferential takes appropriate parameters and returns a function suitable
for use as the third argument to parametric-path. There is some code that
defends against errors (by using ignore-errors) and against certain
peculiarities of IEEE floating-point arithmetic (the code that checks for
not-a-number (NaN) results).

The program is offered here without further comment or apology.

-------------------------------------------------------------------------------

(defparameter units-to-show 4.1)
(defparameter text-width-in-picas 28.0)
(defparameter device-pixels-per-inch 300)
(defparameter pixels-per-unit
  (* (/ (/ text-width-in-picas 6)
        (* units-to-show 2))
     device-pixels-per-inch))

(defparameter big (sqrt (sqrt most-positive-single-float)))
(defparameter tiny (sqrt (sqrt least-positive-single-float)))

(defparameter path-really-losing 1000.0)
(defparameter path-outer-limit (* units-to-show (sqrt 2) 1.1))
(defparameter path-minimal-delta (/ 10 pixels-per-unit))
(defparameter path-outer-delta (* path-outer-limit 0.3))
(defparameter path-relative-closeness 0.00001)
(defparameter back-off-delta 0.0005)

(defun comment-line (stream &rest stuff)
  (format stream "~%% ")
  (apply #'format stream stuff)
  (format t "~%% ")
  (apply #'format t stuff))

(defun parametric-path (from to paramfn plotfn)
  (assert (and (plusp from) (plusp to)))
  (flet ((domainval (x) (funcall paramfn x))
         (rangeval (x) (funcall plotfn (funcall paramfn x)))
         (losing (x) (or (null x)
                         (/= (realpart x) (realpart x))  ;NaN?
                         (/= (imagpart x) (imagpart x))  ;NaN?
                         (> (abs (realpart x)) path-really-losing)
                         (> (abs (imagpart x)) path-really-losing))))
    (when (> to 1000.0)
      (let ((f0 (rangeval from))
            (f1 (rangeval (+ from 1)))
            (f2 (rangeval (+ from (* 2 pi))))
            (f3 (rangeval (+ from 1 (* 2 pi))))
            (f4 (rangeval (+ from (* 4 pi)))))
        (flet ((close (x y)
                 (or (< (careful-abs (- x y)) path-minimal-delta)
                     (< (careful-abs (- x y))
                        (* (+ (careful-abs x) (careful-abs y))
                           path-relative-closeness)))))
          (when (and (close f0 f2)
                     (close f2 f4)
                     (close f1 f3)
                     (or (and (close f0 f1)
                              (close f2 f3))
                         (and (not (close f0 f1))
                              (not (close f2 f3)))))
            (format t "~&Periodicity detected.")
            (setq to (+ from (* (signum (- to from)) 2 pi)))))))
     (let ((fromrange (ignore-errors (rangeval from)))
          (torange (ignore-errors (rangeval to))))
      (if (losing fromrange)
          (if (losing torange)
              '()
              (parametric-path (back-off from to) to paramfn plotfn))
          (if (losing torange)
              (parametric-path from (back-off to from) paramfn plotfn)
              (expand-path (refine-path (list from to) #'rangeval)
                           #'rangeval))))))

(defun back-off (point other)
  (if (or (> point 10.0) (< point 0.1))
      (let ((sp (sqrt point)))
        (if (or (> point sp other) (< point sp other))
            sp
            (* sp (sqrt other))))
      (+ point (* (signum (- other point)) back-off-delta))))

(defun careful-abs (z)
  (cond ((or (> (realpart z) big)
             (< (realpart z) (- big))
             (> (imagpart z) big)
             (< (imagpart z) (- big)))
         big)
        ((complexp z) (abs z))
        ((minusp z) (- z))
        (t z)))

(defparameter max-refinements 5000)

(defun refine-path (original-path rangevalfn)
  (flet ((rangeval (x) (funcall rangevalfn x)))
    (let ((path original-path))
      (do ((j 0 (+ j 1)))
          ((null (rest path)))
        (when (zerop (mod (+ j 1) max-refinements))
              (break "Runaway path"))
        (let* ((from (first path))
               (to (second path))
               (fromrange (rangeval from))
               (torange (rangeval to))
               (dist (careful-abs (- torange fromrange)))
               (mid (* (sqrt from) (sqrt to)))
               (midrange (rangeval mid)))
          (cond ((or (and (far-out fromrange) (far-out torange))
                     (and (< dist path-minimal-delta)
                          (< (abs (- midrange fromrange))
                             path-minimal-delta)
                          ;; Next test is intentionally asymmetric to
                          ;;  avoid problems with periodic functions.
                          (< (abs (- (rangeval (/ (+ to (* from 1.5))
                                                  2.5))
                                     fromrange))
                             path-minimal-delta)))
                 (pop path))
                ((= mid from) (pop path))
                ((= mid to) (pop path))
                (t (setf (rest path) (cons mid (rest path)))))))))
  original-path)

(defun expand-path (path rangevalfn)
  (flet ((rangeval (x) (funcall rangevalfn x)))
    (let ((final-path (list (rangeval (first path)))))
      (do ((p (rest path) (cdr p)))
          ((null p)
           (unless (rest final-path)
             (break "Singleton path"))
           (reverse final-path))
        (let ((v (rangeval (car p))))
          (cond ((and (rest final-path)
                      (not (far-out v))
                      (not (far-out (first final-path)))
                      (between v (first final-path)
                                 (second final-path)))
                 (setf (first final-path) v))
                ((null (rest p))   ;Mustn't omit last point
                 (push v final-path))
                ((< (abs (- v (first final-path))) path-minimal-delta))
                ((far-out v)
                 (unless (and (far-out (first final-path))
                              (< (abs (- v (first final-path)))
                                 path-outer-delta))
                   (push (* 1.01 path-outer-limit (signum v))
                         final-path)))
                (t (push v final-path))))))))

(defun far-out (x)
  (> (careful-abs x) path-outer-limit))

(defparameter between-tolerance 0.000001)

(defun between (p q r)
  (let ((px (realpart p)) (py (imagpart p))
        (qx (realpart q)) (qy (imagpart q))
        (rx (realpart r)) (ry (imagpart r)))
    (and (or (<= px qx rx) (>= px qx rx))
         (or (<= py qy ry) (>= py qy ry))
         (< (abs (- (* (- qx px) (- ry qy))
                    (* (- rx qx) (- qy py))))
            between-tolerance))))

(defun circle (radius)
  #'(lambda (angle) (* radius (cis angle))))

(defun hline (imag)
  #'(lambda (real) (complex real imag)))

(defun vline (real)
  #'(lambda (imag) (complex real imag)))

(defun -hline (imag)
  #'(lambda (real) (complex (- real) imag)))

(defun -vline (real)
  #'(lambda (imag) (complex real (- imag))))

(defun radial (phi quadrant)
  #'(lambda (rho) (repair-quadrant (* rho (cis phi)) quadrant)))

(defun circumferential (rho quadrant)
  #'(lambda (phi) (repair-quadrant (* rho (cis phi)) quadrant)))

;;; Quadrant is 0, 1, 2, or 3, meaning I, II, III, or IV.

(defun repair-quadrant (z quadrant)
  (complex (* (+ (abs (realpart z)) tiny)
              (case quadrant (0 1.0) (1 -1.0) (2 -1.0) (3 1.0)))
           (* (+ (abs (imagpart z)) tiny)
              (case quadrant (0 1.0) (1 1.0) (2 -1.0) (3 -1.0)))))

(defun clamp-real (x)
  (if (far-out x)
      (* (signum x) path-outer-limit)
      (round-real x)))

(defun round-real (x)
  (/ (round (* x 10000.0)) 10000.0))

(defun round-point (z)
  (complex (round-real (realpart z)) (round-real (imagpart z))))

(defparameter hiringshade 0.97)
(defparameter loringshade 0.45)

(defparameter ticklength 0.12)
(defparameter smallticklength 0.09)

;;; This determines the pattern of lines and annuli to be drawn.
(defun moby-grid (&optional (fn 'sqrt) (stream t))
  (comment-line stream "Moby grid for function ~S" fn)
  (shaded-annulus 0.25 0.5 4 hiringshade loringshade fn stream)
  (shaded-annulus 0.75 1.0 8 hiringshade loringshade fn stream)
  (shaded-annulus (/ pi 2) 2.0 16 hiringshade loringshade fn stream)
  (shaded-annulus 3 pi 32 hiringshade loringshade fn stream)
  (moby-lines :horizontal 1.0 fn stream)
  (moby-lines :horizontal -1.0 fn stream)
  (moby-lines :vertical 1.0 fn stream)
  (moby-lines :vertical -1.0 fn stream)
  (let ((tickline 0.015)
        (axisline 0.008))
    (flet ((tick (n) (straight-line (complex n ticklength)
                                    (complex n (- ticklength))
                                    tickline
                                    stream))
           (smalltick (n) (straight-line (complex n smallticklength)
                                         (complex n (- smallticklength))
                                         tickline
                                         stream)))
      (comment-line stream "Real axis")
      (straight-line #c(-5 0) #c(5 0) axisline stream)
      (dotimes (j (floor units-to-show))
        (let ((q (+ j 1))) (tick q) (tick (- q))))
      (dotimes (j (floor units-to-show (/ pi 2)))
        (let ((q (* (/ pi 2) (+ j 1))))
          (smalltick q)
          (smalltick (- q)))))
    (flet ((tick (n) (straight-line (complex ticklength n)
                                    (complex (- ticklength) n)
                                    tickline
                                    stream))
           (smalltick (n) (straight-line (complex smallticklength n)
                                         (complex (- smallticklength) n)
                                         tickline
                                         stream)))
      (comment-line stream "Imaginary axis")
      (straight-line #c(0 -5) #c(0 5) axisline stream)
      (dotimes (j (floor units-to-show))
        (let ((q (+ j 1))) (tick q) (tick (- q))))
      (dotimes (j (floor units-to-show (/ pi 2)))
        (let ((q (* (/ pi 2) (+ j 1))))
          (smalltick q)
          (smalltick (- q)))))))

(defun straight-line (from to wid stream)
  (format stream
          "~%newpath  ~S ~S moveto  ~S ~S lineto  ~S ~
           setlinewidth  1  setlinecap  stroke"
          (realpart from)
          (imagpart from)
          (realpart to)
          (imagpart to)

;;; This function draws the lines for the pattern.
(defun moby-lines (orientation signum plotfn stream)
  (let ((paramfn (ecase orientation
                   (:horizontal (if (< signum 0) #'-hline #'hline))
                   (:vertical (if (< signum 0) #'-vline #'vline)))))
    (flet ((foo (from to other wid)
             (ecase orientation
               (:horizontal
                (comment-line stream
                              "Horizontal line from (~S, ~S) to (~S, ~S)"
                              (round-real (* signum from))
                              (round-real other)
                              (round-real (* signum to))
                              (round-real other)))
               (:vertical
                (comment-line stream
                              "Vertical line from (~S, ~S) to (~S, ~S)"
                              (round-real other)
                              (round-real (* signum from))
                              (round-real other)
                              (round-real (* signum to)))))
             (postscript-path
               stream
               (parametric-path from
                                to
                                (funcall paramfn other)
                                plotfn))
             (postscript-penstroke stream wid)))
      (let* ((thick 0.05)
             (thin 0.02))
        ;; Main axis
        (foo 0.5 tiny 0.0 thick)
        (foo 0.5 1.0 0.0 thick)
        (foo 2.0 1.0 0.0 thick)
        (foo 2.0 big 0.0 thick)
        ;; Parallels at 1 and -1
        (foo 2.0 tiny 1.0 thin)
        (foo 2.0 big 1.0 thin)
        (foo 2.0 tiny -1.0 thin)
        (foo 2.0 big -1.0 thin)
        ;; Parallels at 2, 3, -2, -3
        (foo tiny big 2.0 thin)
        (foo tiny big -2.0 thin)
        (foo tiny big 3.0 thin)
        (foo tiny big -3.0 thin)))))

(defun splice (p q)
  (let ((v (car (last p)))
        (w (first q)))
    (and (far-out v)
         (far-out w)
         (>= (abs (- v w)) path-outer-delta)
         ;; Two far-apart far-out points.  Try to walk around
         ;;  outside the perimeter, in the shorter direction.
         (let* ((pdiff (phase (/ v w)))
                (npoints (floor (abs pdiff) (asin .2)))
                (delta (/ pdiff (+ npoints 1)))
                (incr (cis delta)))
           (do ((j 0 (+ j 1))
                (p (list w "end splice") (cons (* (car p) incr) p)))

;;; This function draws the annuli for the pattern.
(defun shaded-annulus (inner outer sectors firstshade lastshade fn stream)
  (assert (zerop (mod sectors 4)))
  (comment-line stream "Annulus ~S ~S ~S ~S ~S"
                (round-real inner) (round-real outer)
                sectors firstshade lastshade)
  (dotimes (jj sectors)
    (let ((j (- sectors jj 1)))
      (let* ((lophase (+ tiny (* 2 pi (/ j sectors))))
             (hiphase (* 2 pi (/ (+ j 1) sectors)))
             (midphase (/ (+ lophase hiphase) 2.0))
             (midradius (/ (+ inner outer) 2.0))
             (quadrant (floor (* j 4) sectors)))
        (comment-line stream "Sector from ~S to ~S (quadrant ~S)"
                      (round-real lophase)
                      (round-real hiphase)
                      quadrant)
        (let ((p0 (reverse (parametric-path midradius
                                            inner
                                            (radial lophase quadrant)
                                            fn)))
              (p1 (parametric-path midradius
                                   outer
                                   (radial lophase quadrant)
                                   fn))
              (p2 (reverse (parametric-path midphase
                                            lophase
                                            (circumferential outer
                                                             quadrant)
                                            fn)))
              (p3 (parametric-path midphase
                                   hiphase
                                   (circumferential outer quadrant)
                                   fn))
              (p4 (reverse (parametric-path midradius
                                            outer
                                            (radial hiphase quadrant)
                                            fn)))
              (p5 (parametric-path midradius
                                   inner
                                   (radial hiphase quadrant)
                                   fn))
              (p6 (reverse (parametric-path midphase
                                            hiphase
                                            (circumferential inner
                                                             quadrant)
                                            fn)))
              (p7 (parametric-path midphase
                                   lophase
                                   (circumferential inner quadrant)
                                   fn)))
          (postscript-closed-path stream
            (append
              p0 (splice p0 p1) '("middle radial")
              p1 (splice p1 p2) '("end radial")
              p2 (splice p2 p3) '("middle circumferential")
              p3 (splice p3 p4) '("end circumferential")
              p4 (splice p4 p5) '("middle radial")
              p5 (splice p5 p6) '("end radial")
              p6 (splice p6 p7) '("middle circumferential")
              p7 (splice p7 p0) '("end circumferential")
              )))
        (postscript-shade stream
                          (/ (+ (* firstshade (- (- sectors 1) j))
                                (* lastshade j))
                             (- sectors 1)))))))

(defun postscript-penstroke (stream wid)
  (format stream "~%~S setlinewidth   1 setlinecap  stroke"
          wid))

(defun postscript-shade (stream shade)
  (format stream "~%currentgray   ~S setgray   fill   setgray"
          shade))

(defun postscript-closed-path (stream path)
  (unless (every #'far-out (remove-if-not #'numberp path))
    (postscript-raw-path stream path)
    (format stream "~%  closepath")))

(defun postscript-path (stream path)
  (unless (every #'far-out (remove-if-not #'numberp path))
    (postscript-raw-path stream path)))

;;; Print a path as a series of PostScript "lineto" commands.
(defun postscript-raw-path (stream path)
  (format stream "~%newpath")
  (let ((fmt "~%  ~S ~S moveto"))
    (dolist (pt path)
      (cond ((stringp pt)
             (format stream "~%  %~A" pt))
            (t (format stream
                       fmt
                       (clamp-real (realpart pt))
                       (clamp-real (imagpart pt)))
               (setq fmt "~%  ~S ~S lineto"))))))

;;; Definitions of functions to be plotted that are not
;;; standard Common Lisp functions.

(defun one-plus-over-one-minus (x) (/ (+ 1 x) (- 1 x)))

(defun one-minus-over-one-plus (x) (/ (- 1 x) (+ 1 x)))

(defun sqrt-square-minus-one (x) (sqrt (- 1 (* x x))))

(defun sqrt-one-plus-square (x) (sqrt (+ 1 (* x x))))

;;; Because X3J13 voted for a new definition of the atan function,
;;; the following definition was used in place of the atan function
;;; provided by the Common Lisp implementation I was using.

(defun good-atan (x)
  (/ (- (log (+ 1 (* x #c(0 1))))
        (log (- 1 (* x #c(0 1)))))
     #c(0 2)))

;;; Because the first edition had an erroneous definition of atanh,
;;; the following definition was used in place of the atanh function
;;; provided by the Common Lisp implementation I was using.

(defun really-good-atanh (x)
  (/ (- (log (+ 1 x))
        (log (- 1 x)))

;;; This is the main procedure that is intended to be called by a user.
(defun picture (&optional (fn #'sqrt))
  (with-open-file (stream (concatenate 'string
                                       (string-downcase (string fn))
                                       "-plot.ps")
                          :direction :output)
    (format stream "% PostScript file for plot of function ~S~%" fn)
    (format stream "% Plot is to fit in a region ~S inches square~%"
            (/ text-width-in-picas 6.0))
    (format stream
            "%  showing axes extending ~S units from the origin.~%"
            units-to-show)
    (let ((scaling (/ (* text-width-in-picas 12) (* units-to-show 2))))
      (format stream "~%~S ~:*~S scale" scaling))
    (format stream "~%~S ~:*~S translate" units-to-show)
    (format stream "~%newpath")
    (format stream "~%  ~S ~S moveto" (- units-to-show) (- units-to-show))
    (format stream "~%  ~S ~S lineto" units-to-show (- units-to-show))
    (format stream "~%  ~S ~S lineto" units-to-show units-to-show)
    (format stream "~%  ~S ~S lineto" (- units-to-show) units-to-show)
    (format stream "~%  closepath")
    (format stream "~%clip")
    (moby-grid fn stream)
    (format stream
            "~%% End of PostScript file for plot of function ~S"
            fn)
    (terpri stream)))

-------------------------------------------------------------------------------

12.6. Type Conversions and Component Extractions on Numbers

While most arithmetic functions will operate on any kind of number, coercing
types if necessary, the following functions are provided to allow specific
conversions of data types to be forced when desired.

[Function]
float number &optional other

This converts any non-complex number to a floating-point number. With no second
argument, if number is already a floating-point number, then number is
returned; otherwise a single-float is produced. If the argument other is
provided, then it must be a floating-point number, and number is converted to
the same format as other. See also coerce.

[Function]
rational number
rationalize number

Each of these functions converts any non-complex number to a rational number.
If the argument is already rational, it is returned. The two functions differ
in their treatment of floating-point numbers.

rational assumes that the floating-point number is completely accurate and
returns a rational number mathematically equal to the precise value of the
floating-point number.

rationalize assumes that the floating-point number is accurate only to the
precision of the floating-point representation and may return any rational
number for which the floating-point number is the best available approximation
of its format; in doing this it attempts to keep both numerator and denominator
small.

It is always the case that

(float (rational x) x) == x

and

(float (rationalize x) x) == x

That is, rationalizing a floating-point number by either method and then
converting it back to a floating-point number of the same format produces the
original number. What distinguishes the two functions is that rational
typically has a simple, inexpensive implementation, whereas rationalize goes to
more trouble to produce a result that is more pleasant to view and simpler to
compute with for some purposes.

[Function]
numerator rational
denominator rational

These functions take a rational number (an integer or ratio) and return as an
integer the numerator or denominator of the canonical reduced form of the
rational. The numerator of an integer is that integer; the denominator of an
integer is 1. Note that

(gcd (numerator x) (denominator x)) => 1

The denominator will always be a strictly positive integer; the numerator may
be any integer. For example:

(numerator (/ 8 -6)) => -4
(denominator (/ 8 -6)) => 3

There is no fix function in Common Lisp because there are several interesting
ways to convert non-integral values to integers. These are provided by the
functions below, which perform not only type conversion but also some
non-trivial calculations as well.

[Function]
floor number &optional divisor
ceiling number &optional divisor
truncate number &optional divisor
round number &optional divisor

In the simple one-argument case, each of these functions converts its argument
number (which must not be complex) to an integer. If the argument is already an
integer, it is returned directly. If the argument is a ratio or floating-point
number, the functions use different algorithms for the conversion.

floor converts its argument by truncating toward negative infinity; that is,
the result is the largest integer that is not larger than the argument.

ceiling converts its argument by truncating toward positive infinity; that is,
the result is the smallest integer that is not smaller than the argument.

truncate converts its argument by truncating toward zero; that is, the result
is the integer of the same sign as the argument and which has the greatest
integral magnitude not greater than that of the argument.

round converts its argument by rounding to the nearest integer; if number is
exactly halfway between two integers (that is, of the form integer+0.5), then
it is rounded to the one that is even (divisible by 2).

The following table shows what the four functions produce when given various
arguments.

Argument    floor       ceiling     truncate    round
----------------------------------------------------------
 2.6          2           3           2           3
 2.5          2           3           2           2
 2.4          2           3           2           2
 0.7          0           1           0           1
 0.3          0           1           0           0
-0.3         -1           0           0           0
-0.7         -1           0           0          -1
-2.4         -3          -2          -2          -2
-2.5         -3          -2          -2          -2
-2.6         -3          -2          -2          -3
----------------------------------------------------------

If a second argument divisor is supplied, then the result is the appropriate
type of rounding or truncation applied to the result of dividing the number by
the divisor. For example, (floor 5 2) == (floor (/ 5 2)) but is potentially
more efficient.

[change_begin]
This statement is not entirely accurate; one should instead say that (values
(floor 5 2)) == (values (floor (/ 5 2))), because there is a second value to
consider, as discussed below. In other words, the first values returned by the
two forms will be the same, but in general the second values will differ.
Indeed, we have

(floor 5 2) => 2 and 1
(floor (/ 5 2)) => 2 and 1/2

for this example.
[change_end]

The divisor may be any non-complex number.

[change_begin]
It is generally accepted that it is an error for the divisor to be zero.
[change_end]

The one-argument case is exactly like the two-argument case where the second
argument is 1.

[change_begin]
In other words, the one-argument case returns an integer and fractional part
for the number: (truncate 5.3) => 5.0 and 0.3, for example.
[change_end]

Each of the functions actually returns two values, whether given one or two
arguments. The second result is the remainder and may be obtained using
multiple-value-bind and related constructs. If any of these functions is given
two arguments x and y and produces results q and r, then q y+r=x. The first
result q is always an integer. The remainder r is an integer if both arguments
are integers, is rational if both arguments are rational, and is floating-point
if either argument is floating-point. One consequence is that in the
one-argument case the remainder is always a number of the same type as the
argument.

When only one argument is given, the two results are exact; the mathematical
sum of the two results is always equal to the mathematical value of the
argument.

-------------------------------------------------------------------------------
Compatibility note: The names of the functions floor, ceiling, truncate, and
round are more accurate than names like fix that have heretofore been used in
various Lisp systems. The names used here are compatible with standard
mathematical terminology (and with PL/1, as it happens). In Fortran ifix means
truncate. Algol 68 provides round and uses entier to mean floor. In MacLisp,
fix and ifix both mean floor (one is generic, the other flonum-in/fixnum-out).
In Interlisp, fix means truncate. In Lisp Machine Lisp, fix means floor and
fixr means round. Standard Lisp provides a fix function but does not specify
precisely what it does. The existing usage of the name fix is so confused that
it seemed best to avoid it altogether.

The names and definitions given here have recently been adopted by Lisp Machine
Lisp, and MacLisp and NIL (New Implementation of Lisp) seem likely to follow
suit.
-------------------------------------------------------------------------------

[Function]
mod number divisor
rem number divisor

mod performs the operation floor on its two arguments and returns the second
result of floor as its only result. Similarly, rem performs the operation
truncate on its arguments and returns the second result of truncate as its only
result.

mod and rem are therefore the usual modulus and remainder functions when
applied to two integer arguments. In general, however, the arguments may be
integers or floating-point numbers.

(mod 13 4) => 1                 (rem 13 4) => 1
(mod -13 4) => 3                (rem -13 4) => -1
(mod 13 -4) => -3               (rem 13 -4) => 1
(mod -13 -4) => -1              (rem -13 -4) => -1
(mod 13.4 1) => 0.4             (rem 13.4 1) => 0.4
(mod -13.4 1) => 0.6            (rem -13.4 1) => -0.4

-------------------------------------------------------------------------------
Compatibility note: The Interlisp function remainder is essentially equivalent
to the Common Lisp function rem. The MacLisp function remainder is like rem but
accepts only integer arguments.
-------------------------------------------------------------------------------

[Function]
ffloor number &optional divisor
fceiling number &optional divisor
ftruncate number &optional divisor
fround number &optional divisor

These functions are just like floor, ceiling, truncate, and round, except that
the result (the first result of two) is always a floating-point number rather
than an integer. It is roughly as if ffloor gave its arguments to floor, and
then applied float to the first result before passing them both back. In
practice, however, ffloor may be implemented much more efficiently. Similar
remarks apply to the other three functions. If the first argument is a
floating-point number, and the second argument is not a floating-point number
of longer format, then the first result will be a floating-point number of the
same type as the first argument. For example:

(ffloor -4.7) => -5.0 and 0.3
(ffloor 3.5d0) => 3.0d0 and 0.5d0

[Function]
decode-float float
scale-float float integer
float-radix float
float-sign float1 &optional float2
float-digits float
float-precision float
integer-decode-float float

The function decode-float takes a floating-point number and returns three
values.

The first value is a new floating-point number of the same format representing
the significand; the second value is an integer representing the exponent; and
the third value is a floating-point number of the same format indicating the
sign (-1.0 or 1.0). Let b be the radix for the floating-point representation;
then decode-float divides the argument by an integral power of b so as to bring
its value between 1/b (inclusive) and 1 (exclusive) and returns the quotient as
the first value. If the argument is zero, however, the result is equal to the
absolute value of the argument (that is, if there is a negative zero, its
significand is considered to be a positive zero).

The second value of decode-float is the integer exponent e to which b must be
raised to produce the appropriate power for the division. If the argument is
zero, any integer value may be returned, provided that the identity shown below
for scale-float holds.

The third value of decode-float is a floating-point number, of the same format
as the argument, whose absolute value is 1 and whose sign matches that of the
argument.

The function scale-float takes a floating-point number f (not necessarily
between 1/b and 1) and an integer k, and returns (* f (expt (float b f) k)).
(The use of scale-float may be much more efficient than using exponentiation
and multiplication and avoids intermediate overflow and underflow if the final
result is representable.)

Note that

(multiple-value-bind (signif expon sign)
                     (decode-float f)
  (scale-float signif expon))
== (abs f)

and

(multiple-value-bind (signif expon sign)
                     (decode-float f)
  (* (scale-float signif expon) sign))
== f

The function float-radix returns (as an integer) the radix b of the
floating-point argument.

The function float-sign returns a floating-point number z such that z and
float1 have the same sign and also such that z and float2 have the same
absolute value. The argument float2 defaults to the value of (float 1 float1);
(float-sign x) therefore always produces a 1.0 or -1.0 of appropriate format
according to the sign of x. (Note that if an implementation has distinct
representations for negative zero and positive zero, then (float-sign -0.0) =>
-1.0.)

The function float-digits returns, as a non-negative integer, the number of
radix-b digits used in the representation of its argument (including any
implicit digits, such as a ``hidden bit''). The function float-precision
returns, as a non-negative integer, the number of significant radix-b digits
present in the argument; if the argument is (a floating-point) zero, then the
result is (an integer) zero. For normalized floating-point numbers, the results
of float-digits and float-precision will be the same, but the precision will be
less than the number of representation digits for a denormalized or zero
number.

The function integer-decode-float is similar to decode-float but for its first
value returns, as an integer, the significand scaled so as to be an integer.
For an argument f, this integer will be strictly less than

(expt b (float-precision f))

but no less than

(expt b (- (float-precision f) 1))

except that if f is zero, then the integer value will be zero.

The second value bears the same relationship to the first value as for
decode-float:

(multiple-value-bind (signif expon sign)
                     (integer-decode-float f)
  (scale-float (float signif f) expon))
== (abs f)

The third value of integer-decode-float will be 1 or -1.

-------------------------------------------------------------------------------
Rationale: These functions allow the writing of machine-independent, or at
least machine-parameterized, floating-point software of reasonable efficiency.
-------------------------------------------------------------------------------

[Function]
complex realpart &optional imagpart

The arguments must be non-complex numbers; a number is returned that has
realpart as its real part and imagpart as its imaginary part, possibly
converted according to the rule of floating-point contagion (thus both
components will be of the same type). If imagpart is not specified, then
(coerce 0 (type-of realpart)) is effectively used. Note that if both the
realpart and imagpart are rational and the imagpart is zero, then the result is
just the realpart because of the rule of canonical representation for complex
rationals. It follows that the result of complex is not always a complex
number; it may be simply a rational.

[Function]
realpart number
imagpart number

These return the real and imaginary parts of a complex number. If number is a
non-complex number, then realpart returns its argument number and imagpart
returns (* 0 number), which has the effect that the imaginary part of a
rational is 0 and that of a floating-point number is a floating-point zero of
the same format.

[change_begin]
A clever way to multiply a complex number z by i is to write

(complex (- (imagpart z)) (realpart z))

instead of (* z #c(0 1)). This cleverness is not always gratuitous; it may be
of particular importance in the presence of minus zero. For example, if we are
using IEEE standard floating-point arithmetic and z=4+0i, the result of the
clever expression is -0+4i, a true    rotation of z, whereas the result of (* z
#c(0 1)) is likely to be


(4+0i)(+0+i) = ((4)(+0)-(+0)(1))+((4)(1)+(+0)(+0)i
     = ((+0)-(+0))+((4)+(+0))i = +0+4i

which could land on the wrong side of a branch cut, for example.
[change_end]

-------------------------------------------------------------------------------

12.7. Logical Operations on Numbers

The logical operations in this section require integers as arguments; it is an
error to supply a non-integer as an argument. The functions all treat integers
as if they were represented in two's-complement notation.

-------------------------------------------------------------------------------
Implementation note: Internally, of course, an implementation of Common Lisp
may or may not use a two's-complement representation. All that is necessary is
that the logical operations perform calculations so as to give this appearance
to the user.
-------------------------------------------------------------------------------

The logical operations provide a convenient way to represent an infinite vector
of bits. Let such a conceptual vector be indexed by the non-negative integers.
Then bit j is assigned a ``weight''   . Assume that only a finite number of
bits are 1's or only a finite number of bits are 0's. A vector with only a
finite number of one-bits is represented as the sum of the weights of the
one-bits, a positive integer. A vector with only a finite number of zero-bits
is represented as -1 minus the sum of the weights of the zero-bits, a negative
integer.

This method of using integers to represent bit-vectors can in turn be used to
represent sets. Suppose that some (possibly countably infinite) universe of
discourse for sets is mapped into the non-negative integers. Then a set can be
represented as a bit vector; an element is in the set if the bit whose index
corresponds to that element is a one-bit. In this way all finite sets can be
represented (by positive integers), as well as all sets whose complements are
finite (by negative integers). The functions logior, logand, and logxor defined
below then compute the union, intersection, and symmetric difference operations
on sets represented in this way.

[Function]
logior &rest integers

This returns the bit-wise logical inclusive or of its arguments. If no argument
is given, then the result is zero, which is an identity for this operation.

[Function]
logxor &rest integers

This returns the bit-wise logical exclusive or of its arguments. If no argument
is given, then the result is zero, which is an identity for this operation.

[Function]
logand &rest integers

This returns the bit-wise logical and of its arguments. If no argument is
given, then the result is -1, which is an identity for this operation.

[Function]
logeqv &rest integers

This returns the bit-wise logical equivalence (also known as exclusive nor) of
its arguments. If no argument is given, then the result is -1, which is an
identity for this operation.

[Function]
lognand integer1 integer2
lognor integer1 integer2
logandc1 integer1 integer2
logandc2 integer1 integer2
logorc1 integer1 integer2
logorc2 integer1 integer2

These are the other six non-trivial bit-wise logical operations on two
arguments. Because they are not associative, they take exactly two arguments
rather than any non-negative number of arguments.

(lognand n1 n2) == (lognot (logand n1 n2))
(lognor n1 n2) == (lognot (logior n1 n2))
(logandc1 n1 n2) == (logand (lognot n1) n2)
(logandc2 n1 n2) == (logand n1 (lognot n2))
(logorc1 n1 n2) == (logior (lognot n1) n2)
(logorc2 n1 n2) == (logior n1 (lognot n2))

The ten bit-wise logical operations on two integers are summarized in the
following table:

----------------------------------------------------------------
integer1        0       0       1       1
integer2        0       1       0       1       Operation Name
----------------------------------------------------------------
logand          0       0       0       1       and
logior          0       1       1       1       inclusive or
logxor          0       1       1       0       exclusive or
logeqv          1       0       0       1       equivalence (exclusive nor)
lognand         1       1       1       0       not-and
lognor          1       0       0       0       not-or
logandc1        0       1       0       0       and complement of integer1 with integer2
logandc2        0       0       1       0       and integer1 with complement of integer2
logorc1         1       1       0       1       or complement of integer1 with integer2
logorc2         1       0       1       1       or integer1 with complement of integer2

[Function]
boole op integer1 integer2

[Constant]
boole-clr
boole-set
boole-1
boole-2
boole-c1
boole-c2
boole-and
boole-ior
boole-xor
boole-eqv
boole-nand
boole-nor
boole-andc1
boole-andc2
boole-orc1
boole-orc2

The function boole takes an operation op and two integers, and returns an
integer produced by performing the logical operation specified by op on the two
integers. The precise values of the sixteen constants are
implementation-dependent, but they are suitable for use as the first argument
to boole:

----------------------------------------------------------------
integer1        0       0       1       1
integer2        0       1       0       1       Operation Performed
----------------------------------------------------------------
boole-clr       0       0       0       0       always 0
boole-set       1       1       1       1       always 1
boole-1         0       0       1       1       integer1
boole-2         0       1       0       1       integer2
boole-c1        1       1       0       0       complement of integer1
boole-c2        1       0       1       0       complement of integer2
boole-and       0       0       0       1       and
boole-ior       0       1       1       1       inclusive or
boole-xor       0       1       1       0       exclusive or
boole-eqv       1       0       0       1       equivalence (exclusive nor)
boole-nand      1       1       1       0       not-and
boole-nor       1       0       0       0       not-or
boole-andc1     0       1       0       0       and complement of integer1 with integer2
boole-andc2     0       0       1       0       and integer1 with complement of integer2
boole-orc1      1       1       0       1       or complement of integer1 with integer2
boole-orc2      1       0       1       1       or integer1 with complement of integer2

boole can therefore compute all sixteen logical functions on two arguments. In
general,

(boole boole-and x y) == (logand x y)

and the latter is more perspicuous. However, boole is useful when it is
necessary to parameterize a procedure so that it can use one of several logical
operations.

[Function]
lognot integer

This returns the bit-wise logical not of its argument. Every bit of the result
is the complement of the corresponding bit in the argument.

(logbitp j (lognot x)) == (not (logbitp j x))

[Function]
logtest integer1 integer2

logtest is a predicate that is true if any of the bits designated by the 1's in
integer1 are 1's in integer2.

(logtest x y) == (not (zerop (logand x y)))

[Function]
logbitp index integer

logbitp is true if the bit in integer whose index is index (that is, its weight
is   ) is a one-bit; otherwise it is false. For example:

(logbitp 2 6) is true
(logbitp 0 6) is false
(logbitp k n) == (ldb-test (byte 1 k) n)

[change_begin]
X3J13 voted in January 1989 (ARGUMENTS-UNDERSPECIFIED)   to clarify that the
index must be a non-negative integer.
[change_end]

[Function]
ash integer count

This function shifts integer arithmetically left by count bit positions if
count is positive, or right by -count bit positions if count is negative. The
sign of the result is always the same as the sign of integer.

Mathematically speaking, this operation performs the computation   ).

Logically, this moves all of the bits in integer to the left, adding zero-bits
at the bottom, or moves them to the right, discarding bits. (In this context
the question of what gets shifted in on the left is irrelevant; integers,
viewed as strings of bits, are ``half-infinite,'' that is, conceptually extend
infinitely far to the left.) For example:

(logbitp j (ash n k)) == (and (>= j k) (logbitp (- j k) n))

[Function]
logcount integer

The number of bits in integer is determined and returned. If integer is
positive, the 1-bits in its binary representation are counted. If integer is
negative, the 0-bits in its two's-complement binary representation are counted.
The result is always a non-negative integer. For example:

(logcount 13) => 3      ;Binary representation is ...0001101
(logcount -13) => 2     ;Binary representation is ...1110011
(logcount 30) => 4      ;Binary representation is ...0011110
(logcount -30) => 4     ;Binary representation is ...1100010

The following identity always holds:

(logcount x) == (logcount (- (+ x 1)))
             == (logcount (lognot x))

[Function]
integer-length integer

This function performs the computation



This is useful in two different ways. First, if integer is non-negative, then
its value can be represented in unsigned binary form in a field whose width in
bits is no smaller than (integer-length integer). Second, regardless of the
sign of integer, its value can be represented in signed binary two's-complement
form in a field whose width in bits is no smaller than (+ (integer-length
integer) 1). For example:

(integer-length 0) => 0
(integer-length 1) => 1
(integer-length 3) => 2
(integer-length 4) => 3
(integer-length 7) => 3
(integer-length -1) => 0
(integer-length -4) => 2
(integer-length -7) => 3
(integer-length -8) => 3

-------------------------------------------------------------------------------
Compatibility note: This function is similar to the MacLisp function haulong.
One may define haulong as

(haulong x) == (integer-length (abs x))

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

12.8. Byte Manipulation Functions

Several functions are provided for dealing with an arbitrary-width field of
contiguous bits appearing anywhere in an integer. Such a contiguous set of bits
is called a byte. Here the term byte does not imply some fixed number of bits
(such as eight), rather a field of arbitrary and user-specifiable width.

The byte-manipulation functions use objects called byte specifiers to designate
a specific byte position within an integer. The representation of a byte
specifier is implementation-dependent; in particular, it may or may not be a
number. It is sufficient to know that the function byte will construct one, and
that the byte-manipulation functions will accept them. The function byte
accepts two integers representing the position and size of the byte and returns
a byte specifier. Such a specifier designates a byte whose width is size and
whose bits have weights    through   .

[Function]
byte size position

byte takes two integers representing the size and position of a byte and
returns a byte specifier suitable for use as an argument to byte-manipulation
functions.

[Function]
byte-size bytespec
byte-position bytespec

Given a byte specifier, byte-size returns the size specified as an integer;
byte-position similarly returns the position. For example:

(byte-size (byte j k)) == j
(byte-position (byte j k)) == k

[Function]
ldb bytespec integer

bytespec specifies a byte of integer to be extracted. The result is returned as
a non-negative integer. For example:

(logbitp j (ldb (byte s p) n)) == (and (< j s) (logbitp (+ j p) n))

The name of the function ldb means ``load byte.''

-------------------------------------------------------------------------------
Compatibility note: The MacLisp function haipart can be implemented in terms of
ldb as follows:

(defun haipart (integer count)
  (let ((x (abs integer)))
    (if (minusp count)
        (ldb (byte (- count) 0) x)
        (ldb (byte count (max 0 (- (integer-length x) count)))
             x))))

-------------------------------------------------------------------------------

If the argument integer is specified by a form that is a place form acceptable
to setf, then setf may be used with ldb to modify a byte within the integer
that is stored in that place. The effect is to perform a dpb operation and then
store the result back into the place.

[Function]
ldb-test bytespec integer

ldb-test is a predicate that is true if any of the bits designated by the byte
specifier bytespec are 1's in integer; that is, it is true if the designated
field is non-zero.

(ldb-test bytespec n) == (not (zerop (ldb bytespec n)))

[Function]
mask-field bytespec integer

This is similar to ldb; however, the result contains the specified byte of
integer in the position specified by bytespec, rather than in position 0 as
with ldb. The result therefore agrees with integer in the byte specified but
has zero-bits everywhere else. For example:

(ldb bs (mask-field bs n)) == (ldb bs n)

(logbitp j (mask-field (byte s p) n))
   == (and (>= j p) (< j (+ p s)) (logbitp j n))

(mask-field bs n) == (logand n (dpb -1 bs 0))

If the argument integer is specified by a form that is a place form acceptable
to setf, then setf may be used with mask-field to modify a byte within the
integer that is stored in that place. The effect is to perform a deposit-field
operation and then store the result back into the place.

[Function]
dpb newbyte bytespec integer

This returns a number that is the same as integer except in the bits specified
by bytespec. Let s be the size specified by bytespec; then the low s bits of
newbyte appear in the result in the byte specified by bytespec. The integer
newbyte is therefore interpreted as being right-justified, as if it were the
result of ldb. For example:

(logbitp j (dpb m (byte s p) n))
  == (if (and (>= j p) (< j (+ p s)))
         (logbitp (- j p) m)
         (logbitp j n))

The name of the function dpb means ``deposit byte.''

[Function]
deposit-field newbyte bytespec integer

This function is to mask-field as dpb is to ldb. The result is an integer that
contains the bits of newbyte within the byte specified by bytespec, and
elsewhere contains the bits of integer. For example:

(logbitp j (deposit-field m (byte s p) n))
   == (if (and (>= j p) (< j (+ p s)))
          (logbitp j m)
          (logbitp j n))

-------------------------------------------------------------------------------
Implementation note: If the bytespec is a constant, one may of course
construct, at compile time, an equivalent mask m, for example by computing
(deposit-field -1 bytespec 0). Given this mask m, one may then compute

(deposit-field newbyte bytespec integer)

by computing

(logior (logand newbyte m) (logand integer (lognot m)))

where the result of (lognot m) can of course also be computed at compile time.
However, the following expression may also be used and may require fewer
temporary registers in some situations:

(logxor integer (logand m (logxor integer newbyte)))

A related, though possibly less useful, trick is that

(let ((z (logand (logxor x y) m)))
  (setq x (logxor z x))
  (setq y (logxor z y)))

interchanges those bits of x and y for which the mask m is 1, and leaves alone
those bits of x and y for which m is 0.
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

12.9. Random Numbers

The Common Lisp facility for generating pseudo-random numbers has been
carefully defined to make its use reasonably portable. While two
implementations may produce different series of pseudo-random numbers, the
distribution of values should be relatively independent of such
machine-dependent aspects as word size.

[Function]
random number &optional state

(random n) accepts a positive number n and returns a number of the same kind
between zero (inclusive) and n (exclusive). The number n may be an integer or a
floating-point number. An approximately uniform choice distribution is used. If
n is an integer, each of the possible results occurs with (approximate)
probability
1/n. (The qualifier ``approximate'' is used because of implementation
considerations; in practice, the deviation from uniformity should be quite
small.)

The argument state must be an object of type random-state; it defaults to the
value of the variable *random-state*. This object is used to maintain the state
of the pseudo-random-number generator and is altered as a side effect of the
random operation.

-------------------------------------------------------------------------------
Compatibility note: random of zero arguments as defined in MacLisp has been
omitted because its value is too implementation-dependent (limited by fixnum
range).
-------------------------------------------------------------------------------
Implementation note: In general, even if random of zero arguments were defined
as in MacLisp, it is not adequate to define (random n) for integral n to be
simply (mod (random) n); this fails to be uniformly distributed if n is larger
than the largest number produced by random, or even if n merely approaches this
number. This is another reason for omitting random of zero arguments in Common
Lisp. Assuming that the underlying mechanism produces ``random bits'' (possibly
in chunks such as fixnums), the best approach is to produce enough random bits
to construct an integer k some number d of bits larger than (integer-length n)
(see integer-length), and then compute (mod k n). The quantity d should be at
least 7, and preferably 10 or more.

To produce random floating-point numbers in the half-open range [A, B),
accepted practice (as determined by a look through the Collected Algorithms
from the ACM, particularly algorithms 133, 266, 294, and 370) is to compute
X * (B - A) + A, where X is a floating-point number uniformly distributed over
[0.0, 1.0) and computed by calculating a random integer N in the range [0, M)
(typically by a multiplicative-congruential or linear-congruential method mod
M) and then setting X = N/M. See also [27]. If one takes   , where f is the
length of the significand of a floating-point number (and it is in fact common
to choose M to be a power of 2), then this method is equivalent to the
following assembly-language-level procedure. Assume the representation has no
hidden bit. Take a floating-point 0.5, and clobber its entire significand with
random bits. Normalize the result if necessary.

For example, on the DEC PDP-10, assume that accumulator T is completely random
(all 36 bits are random). Then the code sequence

LSH T,-9                 ;Clear high 9 bits; low 27 are random
FSC T,128.               ;Install exponent and normalize

will produce in T a random floating-point number uniformly distributed over
[0.0, 1.0). (Instead of the LSH instruction, one could do

TLZ T,777000             ;That's 777000 octal

but if the 36 random bits came from a congruential random-number generator, the
high-order bits tend to be ``more random'' than the low-order ones, and so the
LSH would be better for uniform distribution. Ideally all the bits would be the
result of high-quality randomness.)

With a hidden-bit representation, normalization is not a problem, but dealing
with the hidden bit is. The method can be adapted as follows. Take a
floating-point 1.0 and clobber the explicit significand bits with random bits;
this produces a random floating-point number in the range [1.0, 2.0). Then
simply subtract 1.0. In effect, we let the hidden bit creep in and then
subtract it away again.

For example, on the DEC VAX, assume that register T is completely random (but a
little less random than on the PDP-10, as it has only 32 random bits). Then the
code sequence

INSV #^X81,#7,#9,T     ;Install correct sign bit and exponent
SUBF #^F1.0,T          ;Subtract 1.0

will produce in T a random floating-point number uniformly distributed over
[0.0, 1.0). Again, if the low-order bits are not random enough, then the
instruction

ROTL #7,T

should be performed first.

Implementors may wish to consult reference [41] for a discussion of some
efficient methods of generating pseudo-random numbers.
-------------------------------------------------------------------------------

[Variable]
*random-state*

This variable holds a data structure, an object of type random-state, that
encodes the internal state of the random-number generator that random uses by
default. The nature of this data structure is implementation-dependent. It may
be printed out and successfully read back in, but may or may not function
correctly as a random-number state object in another implementation. A call to
random will perform a side effect on this data structure. Lambda-binding this
variable to a different random-number state object will correctly save and
restore the old state object.

[Function]
make-random-state &optional state

This function returns a new object of type random-state, suitable for use as
the value of the variable *random-state*. If state is nil or omitted,
make-random-state returns a copy of the current random-number state object (the
value of the variable *random-state*). If state is a state object, a copy of
that state object is returned. If state is t, then a new state object is
returned that has been ``randomly'' initialized by some means (such as by a
time-of-day clock).

-------------------------------------------------------------------------------
Rationale: Common Lisp purposely provides no way to initialize a random-state
object from a user-specified ``seed.'' The reason for this is that the number
of bits of state information in a random-state object may vary widely from one
implementation to another, and there is no simple way to guarantee that any
user-specified seed value will be ``random enough.'' Instead, the
initialization of random-state objects is left to the implementor in the case
where the argument t is given to make-random-state.

To handle the common situation of executing the same program many times in a
reproducible manner, where that program uses random, the following procedure
may be used:

  1.  Evaluate (make-random-state t) to create a random-state object.

  2.  Write that object to a file, using print, for later use.

  3.  Whenever the program is to be run, first use read to create a copy of the
     random-state object from the printed representation in the file. Then use
     the random-state object newly created by the read operation to initialize
     the random-number generator for the program.

It is for the sake of this procedure for reproducible execution that
implementations are required to provide a read/print syntax for objects of type
random-state.

It is also possible to make copies of a random-state object directly without
going through the print/read process, simply by using the make-random-state
function to copy the object; this allows the same sequence of random numbers to
be generated many times within a single program.
-------------------------------------------------------------------------------
Implementation note: A recommended way to implement the type random-state is
effectively to use the machinery for defstruct. The usual structure syntax may
then be used for printing random-state objects; one might look something like

#S(RANDOM-STATE DATA #(14 49 98436589 786345 8734658324 ...))

where the components are of course completely implementation-dependent.
-------------------------------------------------------------------------------

[Function]
random-state-p object

random-state-p is true if its argument is a random-state object, and otherwise
is false.

(random-state-p x) == (typep x 'random-state)

-------------------------------------------------------------------------------

12.10 Implementation Parameters

The values of the named constants defined in this section are
implementation-dependent. They may be useful for parameterizing code in some
situations.

[Constant]
most-positive-fixnum
most-negative-fixnum

The value of most-positive-fixnum is that fixnum closest in value to positive
infinity provided by the implementation.

The value of most-negative-fixnum is that fixnum closest in value to negative
infinity provided by the implementation.

[change_begin]
X3J13 voted in January 1989 (FIXNUM-NON-PORTABLE)   to specify that fixnum must
be a supertype of the type (signed-byte 16), and additionally that the value of
array-dimension-limit must be a fixnum. This implies that the value of
most-negative-fixnum must be less than or equal to   , and the value of
most-positive-fixnum must be greater than or equal to both    and the value of
array-dimension-limit.
[change_end]

[Constant]
most-positive-short-float
least-positive-short-float
least-negative-short-float
most-negative-short-float

The value of most-positive-short-float is that short-format floating-point
number closest in value to (but not equal to) positive infinity provided by the
implementation.

The value of least-positive-short-float is that positive short-format
floating-point number closest in value to (but not equal to) zero provided by
the implementation.

The value of least-negative-short-float is that negative short-format
floating-point number closest in value to (but not equal to) zero provided by
the implementation. (Note that even if an implementation supports minus zero as
a distinct short floating-point value, least-negative-short-float must not be
minus zero.)

[change_begin]
X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to clarify that these definitions
are to be taken quite literally. In implementations that support denormalized
numbers, the values of least-positive-short-float and
least-negative-short-float may be denormalized.
[change_end]

The value of most-negative-short-float is that short-format floating-point
number closest in value to (but not equal to) negative infinity provided by the
implementation.

[Constant]
most-positive-single-float
least-positive-single-float
least-negative-single-float
most-negative-single-float
most-positive-double-float
least-positive-double-float
least-negative-double-float
most-negative-double-float
most-positive-long-float
least-positive-long-float
least-negative-long-float
most-negative-long-float

These are analogous to the constants defined above for short-format
floating-point numbers.

[change_begin]
[Constant]
least-positive-normalized-short-float
least-negative-normalized-short-float

X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to add these constants to the
language.

The value of least-positive-normalized-short-float is that positive normalized
short-format floating-point number closest in value to (but not equal to) zero
provided by the implementation. In implementations that do not support
denormalized numbers this may be the same as the value of
least-positive-short-float.

The value of least-negative-normalized-short-float is that negative normalized
short-format floating-point number closest in value to (but not equal to) zero
provided by the implementation. (Note that even if an implementation supports
minus zero as a distinct short floating-point value,
least-negative-normalized-short-float must not be minus zero.) In
implementations that do not support denormalized numbers this may be the same
as the value of least-positive-short-float.

[Constant]
least-positive-normalized-single-float
least-negative-normalized-single-float
least-positive-normalized-double-float
least-negative-normalized-double-float
least-positive-normalized-long-float
least-negative-normalized-long-float

These are analogous to the constants defined above for short-format
floating-point numbers.
[change_end]

[Constant]
short-float-epsilon
single-float-epsilon
double-float-epsilon
long-float-epsilon

These constants have as value, for each floating-point format, the smallest
positive floating-point number e of that format such that the expression

(not (= (float 1 e) (+ (float 1 e) e)))

is true when actually evaluated.

[Constant]
short-float-negative-epsilon
single-float-negative-epsilon
double-float-negative-epsilon
long-float-negative-epsilon

These constants have as value, for each floating-point format, the smallest
positive floating-point number e of that format such that the expression

(not (= (float 1 e) (- (float 1 e) e)))

is true when actually evaluated.

-------------------------------------------------------------------------------